{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1-异常检测"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# note:\n",
    "* [covariance matrix](http://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html)\n",
    "* [multivariate_normal](http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.multivariate_normal.html)\n",
    "* [seaborn  bivariate kernel density estimate](https://stanford.edu/~mwaskom/software/seaborn/generated/seaborn.kdeplot.html#seaborn.kdeplot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\sklearn\\cross_validation.py:41: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(context=\"notebook\", style=\"white\", palette=sns.color_palette(\"RdBu\"))\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.io as sio\n",
    "from scipy import stats\n",
    "from sklearn.cross_validation import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You want to divide data into 3 set. \n",
    "1. Training set\n",
    "2. Cross Validation set\n",
    "3. Test set.  \n",
    "\n",
    "You shouldn't be doing prediction using training data or Validation data as it does in the exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['__header__', '__version__', '__globals__', 'X', 'Xval', 'yval'])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mat = sio.loadmat('./data/ex8data1.mat')\n",
    "mat.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = mat.get('X')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "divide original validation data into validation and test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Xval, Xtest, yval, ytest = train_test_split(mat.get('Xval'),\n",
    "                                            mat.get('yval').ravel(),\n",
    "                                            test_size=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFXCAYAAAC7nNf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPW97//XmmvITJJJzA0IIKjU24Zw89T9OOKjWopH\nZSuKRaqwz4bdh8djVWwrFVpbWhSs1Xbvw/G+fTzsRru9UNtStWhBe9w/dXPThCKi3C0JuQCZJDPJ\n3Nfvj5UMmRBiQCYhi/fzL2Zmzcx3ZUje8/2u7/fzNUzTNBEREZFBzzHQDRAREZFTQ6EuIiJiEwp1\nERERm1Coi4iI2IRCXURExCZcA92ALyMSibBt2zZKSkpwOp0D3RwREZGsSiaTNDY2cvHFF5OTk3PM\n44M61Ldt28Ytt9wy0M0QERHpVy+88AKTJ08+5v5BHeolJSWAdXLl5eUD3BoREZHsqqur45Zbbknn\nX3dZC/V4PM6SJUuoqakhFotx++23M3ToUG677TbOPvtsAObMmcPVV1+dfk4qlWLp0qV8+umneDwe\nHnjgAUaNGnXc9+gcci8vL6eioiJbpyIiInJaOd4l56yF+po1awgEAvziF78gGAxy/fXXc8cdd/BP\n//RPzJ8/v8fnrFu3jlgsxksvvURVVRUPPfQQTzzxRLaaKCIiYitZC/WrrrqK6dOnA2CaJk6nk23b\ntrF3717Wr1/PqFGjWLJkCX6/P/2cLVu2cNlllwFQWVnJtm3bstU8ERER28nakjafz4ff7ycUCnHX\nXXexcOFCxo0bx6JFi3jhhRcYMWIEjz32WMZzQqFQRsg7nU4SiUS2migiImIrWV2nfvDgQebNm8d1\n113HjBkzmDZtGhdffDEA06ZNY/v27RnH+/1+wuFw+nYqlcLlGtRz+URERPpN1kL90KFDzJ8/n3vv\nvZdZs2YBsGDBArZu3QrABx98wEUXXZTxnIkTJ/Luu+8CUFVVxdixY7PVPBEREdvJWjf4ySefpKWl\nhccff5zHH38cgPvuu4/ly5fjdrspLi5m2bJlACxatIiFCxcybdo03nvvPW6++WZM02T58uXZap6I\niIjtGIN5P/UDBw5w5ZVXsn79ei1pExER2/ui3FPtdxEREZtQqIuIiNiEppZ3MFMpWqq3Em1owFta\nSn7leAzDGOhmiYiI9JlCvUNL9VaaNm8GoO3zzwEomFA5kE0SERE5IRp+7xBtaOj1toiIyOlOod7B\nW1ra620REZHTnYbfO+RXjgfIuKYuIiIymCjUOxiGoWvoIiIyqGn4XURExCYU6iIiIjahUBcREbEJ\nhbqIiIhNKNRFRERsQqEuIiJiEwp1ERERm1Coi4iI2IRCXURExCYU6iIiIjahUBcREbEJhbqIiIhN\nKNRFRERsQqEuIiJiEwp1ERERm1Coi4iI2IRCXURExCZc2XrheDzOkiVLqKmpIRaLcfvttzNs2DCW\nLVuG0+nE4/Hw85//nOLi4oznzZw5E7/fD0BFRQUrVqzIVhNFRERsJWuhvmbNGgKBAL/4xS8IBoNc\nf/31VFRUcP/993PBBRfw4osv8swzz7B48eL0c6LRKKZpsmrVqmw1S0RExLayFupXXXUV06dPB8A0\nTZxOJ7/85S8pLS0FIJlM4vV6M56zY8cO2tvbmT9/PolEgu9+97tUVlZmq4kiIiK2krVQ9/l8AIRC\nIe666y4WLlyYDvQPP/yQ559/nhdeeCHjOTk5OSxYsICbbrqJffv28e1vf5u1a9ficmWtmSIiIraR\n1bQ8ePAgd9xxB9/61reYMWMGAG+88QZPPPEETz/9NEVFRRnHjx49mlGjRmEYBqNHjyYQCNDY2MjQ\noUOz2UwRERFbyNrs90OHDjF//nzuvfdeZs2aBcAf/vAHnn/+eVatWsWIESOOec7q1at56KGHAKiv\nrycUClFSUpKtJoqIiNhK1nrqTz75JC0tLTz++OM8/vjjJJNJdu7cybBhw7jzzjsBmDJlCnfddReL\nFi1i4cKFzJo1i8WLFzNnzhwMw2D58uUaehcREekjwzRNc6AbcbIOHDjAlVdeyfr166moqBjo5oiI\niGTVF+Weis+IiIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGxCoS4iImITCnURERGb\nUKiLiIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGxCoS4iImITCnURERGbUKiLiIjY\nhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGxCoS4iImITCnURERGbUKiLiIjYhEJdRETE\nJlzZeuF4PM6SJUuoqakhFotx++23c+6553LfffdhGAbnnXceP/nJT3A4jn6vSKVSLF26lE8//RSP\nx8MDDzzAqFGjstVEERERW8laT33NmjUEAgF+85vf8G//9m8sW7aMFStWsHDhQn7zm99gmibr16/P\neM66deuIxWK89NJLfO973+Ohhx7KVvNERERsJ2uhftVVV3H33XcDYJomTqeTjz/+mEsuuQSAqVOn\n8v7772c8Z8uWLVx22WUAVFZWsm3btmw1T0RExHayFuo+nw+/308oFOKuu+5i4cKFmKaJYRjpx1tb\nWzOeEwqF8Pv96dtOp5NEIpGtJoqIiNhKVifKHTx4kHnz5nHdddcxY8aMjOvn4XCY/Pz8jOP9fj/h\ncDh9O5VK4XJl7bK/iIiIrWQt1A8dOsT8+fO59957mTVrFgAXXnghGzZsAODdd99l8uTJGc+ZOHEi\n7777LgBVVVWMHTs2W80TERGxnayF+pNPPklLSwuPP/44c+fOZe7cuSxcuJCVK1cye/Zs4vE406dP\nB2DRokXU1tYybdo0PB4PN998MytWrGDx4sXZap6IiIjtGKZpmgPdiJN14MABrrzyStavX09FRcVA\nN0dERCSrvij3VHxGRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGxCoS4iImITCnURERGbUKiL\niIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGxCoS4iImITCnURERGbUKiLiIjYhEJd\nRETEJhTqIiIiNqFQFxERsQmFuoiIiE24BroBIiJyfGYqRUv1VqINDXhLS8mvHI9hGAPdLDlNKdRF\nRE5jLdVbadq8GYC2zz8HoGBC5UA2SU5jCnUROWnqRWZftKGh19siXSnUReSkqReZfd7S0vTPtvO2\nyPEo1EXkpKkXmX35leMBMkZDRI4nq6FeXV3NI488wqpVq7jnnns4dOgQADU1NYwfP55f/epXGcfP\nnDkTv98PQEVFBStWrMhm80TkS1IvMvsMw9Doh/RZ1kL9mWeeYc2aNQwZMgQgHeDNzc3MmzePxYsX\nZxwfjUYxTZNVq1Zlq0kicoqpFylyesnaOvWRI0eycuXKY+5fuXIlt956K6XdvtHv2LGD9vZ25s+f\nz7x586iqqspW00TkFOnsRZZO/wYFEyo1SU5kgGUt1KdPn47LlTkQcPjwYT744ANuuOGGY47Pyclh\nwYIFPPvss/z0pz/l+9//PolEIlvNExERsZ1+nSi3du1arr32WpxO5zGPjR49mlGjRmEYBqNHjyYQ\nCNDY2MjQoUP7s4kiIiKDVr+Wif3ggw+YOnVqj4+tXr2ahx56CID6+npCoRAlJSX92TwR6SMzlaL5\noyoa3nyL5o+qME1zoJskIvRzqO/du5cRI0Zk3Ldo0SJqa2uZNWsWra2tzJkzh3vuuYfly5cfM3wv\nIqeHzvXpbZ9/TtPmzbRUVZ+y19YXBpGTl9XUrKio4OWXX07ffv3114855uGHH07/+9FHH81mc0Tk\nFMnm+nQVtBE5edqlTUROWPf16KdyfboK2oicPI1vi5wBTnWN9myuT1dBG5GTp1AXOQOc6iHtbFY5\nU0EbkZOnUBcZhE60532qhrS/6H1PxYjAF31h0M5wIsenUBcZhE60532qhrSbP6qi7o03SIbCOP0+\nTNMkMHHCSbfreHoLbk2kEzk+hbrIIHSiPe9TNaQd3LyFSE0tAPHmZoKbNmeE+qkaEegtuDWRTuT4\nFOoig9CJ9ryzdQ08UldHw5tvpb8onKoRga5BbZomTRs3pb+QeEpLNJFO5DgU6iKDUH9MJutpCDxw\nyWSi9fUkwmFS0SgYBm2ff54O2d7adbwh9Z7u7/rlIFpfj2ECDuu9ApMmUTh5MtGGBjwlJaSSSfY9\n9QwAgUsmU1CpjWXkzKVQFxmEstXzNlMpmquqCW7aTKSuDoCcoUOPDoFXVmJgEG1ooG3ffswu2Rlt\naOi1XccbUu/p/q5fDkiZ0CWkY42NlEz7uvW8jZto+eQTku3tGIZhfQFA+4/LmUuhLnIG+aKZ4y3V\nW6l/7XXaa2uJH2nC8LgxDAPv0PJjQrv5o6p0GEPPw+CdJV+Dm7cQ3r0bp9+Ht7wcwzCsLw0fVdGw\nbh3JSARvWZkVzH14n84vAuE9e4jWHsTwuHH5/STC4fTr9nSOmjkvdqdQFzmDfNHM8WhDA4lwGADD\n48aMxUmEw3g5NrS7D7Xnjfs7mruFaUv1VureeIP2AzXEDx0iFU+QO2Y0eRdeQKI1RNPfNhGpqSVS\ne5DY4SPkXXjBcd8nUldHojVEpL6e9v2fg2Hg8vnS7QRw+XzW6x7oOMf9+wnt3IW7IB9vaSkmJsHN\nW457/iKn0kB8iVSoi9hY9z8qkfr6jMe7zxz3lpbi8vmINzfj8vtx5Hhx5uZCyiSVTBL86CNiDY3p\nP1BdA7Frjzp9PbyhgWQoTDIcJpVIAJAIBvEECnHl59G6YwfRQ4dIRiJEaw9y1t9fesz8gHSv/aMq\nmg5sJt4cJNJgXWf3lpdjmiaxI0dwer34z/8Kzjw/8eag9f519YR378F37jm07d9P5GAdybY2XD5r\nxEAz5yWbBmL5pUJdxMa6/1FxFwSOPmiaxJtbMmav51eOt3qzG63nuAIBKyANg/o//ckK0qHlGX+g\nOr84NKxbR6o9gre8HIyjPXin35fuSbvz8/CNGZPuOUdqakm0toJhjQzEm4LH9GR6en1vWRmGCbmj\nRuIOBIgFmzAMg3hzM3D0+YmwtZ4erICP1NZC+jgonDI5Kz93ERiY5Zd9CvWdO3dy3nnnZdxXVVVF\nZaWGrUROZ93/iLjy8/Cfdy7RhgbizS3Egk3Em4MZIR2YMIHABGvtecObb5FosQIwGbKG5b1kLjPr\nfJ1kJEKktpbYkSM4vF7cBQGKp33d2jo1ZRKpPUjO8GFWKHd8icgZNpR4k3Xt3unz9XgOnQVvIjU1\nJEJhTNMkZ9hQCqdMpmBCJfV/Wktoxw4SoRCpWAzv0KH4Ro/GlZ+Hu8AKfLAC3jt8GIZhkAyF8ZaV\nqQStZNVA7GPQa6hv2bKFVCrFj370Ix588MH0vsaJRIKlS5fy5ptvZr2BInLyuv9RySkrSw//Nbz5\nVnqYGnruRXR9vtPvs5aWkbnMLLyrYwJcWRmxw0eIB4PknjOGWLCJ1uqtBCZOsGa5V1VnXH9vqaq2\nntPUhMPjwZXnJ9BDz/lowRurB54MhymcPDkdyPGWVlo/+4xkSyvJeIyUaeIpKsR/3rnkT/9G+n07\nA75zJKBw8mRNkpOsGoh9DHoN9ffff5+NGzfS0NDAv/7rvx59ksvF7Nmzs944EflyevujcrxeRNfr\n8J7SEgomTKB5y4d4S0rxFBXiys+HlIkJRA/WEWsKkjxwgGQoTKq9ndzRo8kpL0+/Lxy7BK/z+rvh\ncOAtKsJbVkbhJVOs4f9u8wA6OxMALr8f3znnZLxW7MgRUu3t1gQ/wyDVHkm/d9f3NU0z44uFeumS\nbdnc+Oh4eg31O++8E4Df//73XH/99f3SIBHpm77MrDUMg/zx49LHtVRVp487XuD3eB3eYWBgEGtq\nov1ADdG6OqKHDmG43cQOHSLZ1gYYeEuKM9aUe0tLe2xnelTAsK7R544ceXQJW1V1xvt7igrJGTqU\nSO1BANyFAUzTTJ9rrKGBVCKB4XGTisZIdcze7z7UORB/YEX6W5+uqW/YsIENGzYcc/+KFStOeYNE\npG/6OrP2eMd1DblUIkHtK7+lbf8+km3t1lpyhxWabfv34Q5YE+xCn+wgUnsQV2GA2KHDGC4XZiKB\nw+PBjMdxFxWRU1ZG7siR6WH22ldW0/jOO2CCd9hQWj/bSSIYJFpfn55U1zWAu18GcPr9OH1+Iget\nUD/03vvknnsOhRMnAuAtL8OVl0cqFsPMTeLwWbP1TdPMCH+RM0GfQv2SSy5J/zuRSLB+/XrGjBmT\ntUaJSO/MVIqmjZsI79mDy+fDU1aWUR+9a6+9tzrqnccdfPX3NL79NgCJUIh4UxN5F10IQO6os9PX\n3uPBoNUj7ugNp6IRXIEAydYQqWSCtj17yDv/fEq+MQ1Mk9rVr1L721eJt7bi8HiIHjpEtPYgBZMm\nWtXoTJPCKVN6vSyQDIVp/uhDEqEQAOHPPqP+j6+nQ71wypSjS+eiUTyFheAwCG7ZAqaJ4XCo2Iyc\nMfoU6jNnzsy4PWvWLObMmZOVBolI3yq/RRrqiTc3E29uJnbkCO6zitL10eFor/14ddTD+/bR+M5f\nSCXitFRvJRWP4/R4cPn8ODweckeOxF1URGj3boKbrIIt7rIyIgcOWAFrWNe4XT4fTq+XZCicXl7W\nUlUNQNOWzSRaWkm2tWHGEzi8HsAaCveWlYFpFZUJvbILV34eOWVl5I0fZ7W1y9p6MxaHZBLTNElB\nukBO53l2VqJr27f/6PC/aVL32usZ69K7/lxE7Oik1qnv3r2bBhVtEMma7kPmZiqV0eOM1NdboUjH\nUjPTTN+GzN758eqodw6lGx430bp6zFQSh9uDK6+dkiu/Run0b1Dz8moaXnudeGsrAEOG5OAtLibh\nduMKBPCffz4OwyDe0kyyrd3a4GX3HgwThowaSSoawzSsEDcTCVxnFZEzfLjVno4vGJGGeiI1tThz\ncnB4vQR27qJ85nWEfreL5q1bMZwujNxcrOlyVs/bW1p6TPU6wzAyCuBEDtYR2vkZZsy63m6aporN\niO31KdTPP/9865eyYxZqUVER3/3ud7PaMJEzWffwCW7eAp3XuDsmrxmGkZ5l7i4IZCxPy5gk1mX2\nuDsQoPXTT0mGw1YhFrdVYtV0OTFjVk84EQ7TvG07mK9w+D//P+Itrel6LsnmZoouvTTdFoDApEmE\nd+2mcf36jh68QTwYJBEKEWtttV7fNHH5fASmXkZucQltn++HlImnvJy2PXusIf8jTbiLCglu2UL7\ngQOEPvvsaLsL8mHEcFLROE6/j0hdPQdfew1vWVnGyETXLzCR2oOYJsRbWzCTSTBNymdcewo+HZHT\nV59CfceOHdluh4h00f26cnddi8h4S0vJGz+O1m7D9Z269vrba2pp//xzSKWg44u6w+PGYYLD5yOV\nTFrPqaqibc/ujglwMUyw1pIHAgSmTE5vyBJvbmH/c78mevAgqWQSM57A6XGD30/os89IHDpEKh7v\n6BSkiB2owe3x4C4o6Kj5XotzyBDMWByHxw1Y6+Hb9u7NPN+cHPKnTqV1+ydEaw8SjR0k0doCQE6X\ncq9dJ/+17dtPaOdO634MrLJ1X/6zETmd9SnUW1paWLlyJf/1X/+Fy+Vi6tSp3H777eTk5GS7fSJn\npO7LzbpuRAKZRWQ6He9acWfgmaZJeOcu4sEg7sIArkAAZ04OOcOHEW8KEjl4kPgRq/qaAaRiMRxD\ncnEVFpEMtZJTVk7hV79KfuV4HA4H5pYPqXllNW179lozzTExMEh5PDiGDMGZm4uZSGC4rD8zDreH\nZKs12S1aV08yEukoD+sh95wxkEpZRWxKS0m2tRPdvSe9+1pg8iTAIBEM4vC4MbGK0KSr3PVQqStw\nyWSCH36Y/sLgGzuWWEPjSX8mIoNBn0L93nvvZcyYMTzyyCOYpslvf/tbfvjDH/Loo49mu30iZ6Tu\na6pN00zvY+4pKcE0zYya7b3N6O7s9Ufr60m0tmKmkla9daBw4kTOvu3bmKZJ7cuvcHDNa9Zjbhdm\nMkmqrQ0zGgEMnLlDSLa00Fq9lYIJlQQ3byFaXw+miRm3arubhoHpdBAPBkm2tZOKJwATMHC4XeSO\nGQ2mSXtNLfGmJrzDhuI79xxyR44kp6wsXXbWNCHV3k48GMQ/dixDZ93IoT+vwzdmDO21tZ0/FHxj\nxmRUl+uqoLKSkiuvJLhlS3qiXG9lOrUtq9hBn0K9pqaGp556Kn37hz/8Iddeq2tTItkOgu6vX/KN\nabR0K84Cx/bSM6rClZQQmDyJxnXr8Y0ZTbTxEKlYzApjTJo/qiK/cjy+887D/5WxRA7UEO9YPmYA\nqXgCE6t2u6eoiGhDA2YqRaSuDjOZsiawGQaYJk5fLobDiRmLYeTm4i4+CzMSweFykXfRRYy++07q\nfvs7wMTwuEm2txOtr6doypR0kZzmqmraDx60ZuCfM4a8Cy/A0TE5rnMGeyIcJjBpEsNuurHHn3fn\n+bvy8whMmpSeWd9bFbmB2FFL5FTrU6iPGjWKzZs3M3myVZd5x44djBo16gufV11dzSOPPMKqVavY\nvn07t912G2effTYAc+bM4eqrr04fm0qlWLp0KZ9++ikej4cHHnigT+8hMpCyHQQ9vX7X4fRofT0N\n69YBHUP2pmk9Z+OmdHGXts8/p3DyZEq//nWaNm3GU1dHe00NziFDAIO6P75G08ZNgFXIxQCorQVM\n4keCONyudMW2RDiMt7SUluqtYBh4zjqL6KFDpKJRq2cfjYLDidPrtcI2kbS+VEyaROGUyTidTtwF\n+QQmTSJaV0ciHCan1Arbzi8rkYMHadu7F1deHi6/n8RXWo+eH/TpC1TXnxuA/7xzv/BzGYgdtURO\ntT6F+ueff86tt97K6NGjcTqd7N27l4KCAq644goMw2D9+vXHPOeZZ55hzZo1DBkyBICPP/6Yf/qn\nf2L+/Pk9vse6deuIxWK89NJLVFVV8dBDD/HEE098iVMTyb5sB0FPr991OD1SU8uQYcMyAqxp82bC\ne/aktxf1DrUmkpV8YxqmaRLctJlEKITT77deo7aWRDiMy+fDNKxr1QDe8qHEm4LEW1oxnA5isTgO\nt5tkKkW8sRFveZnVlf/UJHboiDW0H0+Cw4HpcJBoa7dmrQPO/HxCO3fSsG4dhtOFw5eLd2g5Xo5u\nrNJ5rg6vF3deHobLxZBhw3Dl5wEnVub1ZD6XgdhRS+RU61OoP/nkkyf8wiNHjmTlypUsWrQIgG3b\ntrF3717Wr1/PqFGjWLJkCX6/P338li1buOyyywCorKxk27ZtJ/yeIv0t20HQ0+t39lgb3vozzpwc\nqxDLwToi9fXpnqvT7yPe3Gz1rDueZxiG9bjDwJnnt3Y+sy534/T78JSWET14kHhTkERbiOS+KIn2\n9vRWpYbLRezwEfY99ji5o0ZhYhI/0kSs8bBV+z2VshqZTJJsb8c0TTzFZ5E/cQLhTz9Nl3kF8I8d\ni/+iC9PnY6ZSxJtbCO/aTSoaxenz4xhinVuipfWEy72ezOcyEDtqiZxqfQr10tJS3n//fZqamjLu\n722Tl+nTp3PgwIH07XHjxnHTTTdx8cUX88QTT/DYY4/xgx/8IP14KBTKCHmn00kikcDlOqn6OCL9\nIttBcLwh54IJlYR27iK8dy/JSIR4czP+llb8551L235rDTimiZGTQ6K1lYY/ryO0cxfOPOt3LF24\npjWEM8+Pt6yM6ME6WrZvJx5shlQSDIdVU92dTK9LT8Vj1rB+SQnRxkar/rvbqv/ek9ihw4Q+2UGi\nKUgqmcDp81nL25IJSqd/I31cc1U1sWATTr8PExP3kCGYySQun494MEhLVfUJXdY4mc9FG76IHfQp\nMe+++24aGxs555xzMr4tn8jObdOmTSM/Pz/972XLlmU87vf7CXcp/ZhKpRToctrLZhB80SQ8V34e\nOcOHkQyFcfp9uPLzyK8cT2jnLkK7duMYMoTwrl2EIhE8JSWE9+zBP3Ysrjx/unBN4JpJ6aHv1o+3\nE2u01pUDODo3a8nJgVAYnE6rR55M0VxdjeH24BySg7uwkFQsTrKtDcPhwEwkcObm4nC7MAxItUdw\nFwYId6w9d/n95I46O+NcO7dJzSkvxzRN2nbttsLfAJPjV4I73s9IAS1nqj6l5p49e1i7du2XeqMF\nCxZw//33M27cOD744AMuuuiijMcnTpzIO++8w9VXX01VVRVjx479Uu8nMpiZqRS1q18luGWLtXa7\no2fdNahyysrSFeU6bxuGgbsgH1een/YDNcRqD5KKx0lFo3iHD8dMJCicPNmaFV9cDMbR3qzD682o\nm24ChtMJDgeeslJSsTjxw4chlSLRGsJwOMDMx1NURO7oszEcDpJtbcSDzTi8HqtYjc/HkOHD8JRb\n7TeTSc667DKGzrohHciRujpaP9lBpK4OV54fTJNULJYegQAomjKlx5/TF01U1DI1OdP0KdRHjhxJ\nbW0tw4YNO+k3Wrp0KcuWLcPtdlNcXJzuqS9atIiFCxcybdo03nvvPW6++WZM02T58uUn/V4ig0Fv\ngdNSvZXgli3pDVvih4+Q6Ai4L9oP3VNcTKSmhvb9n5PsWLoWCzaTSiTxFBURqa8np6yMZCLOgVX/\nQSIYxBUIkDO8AndhgHhTEDOZxOF2k4rHSIVjuPLycLrdpHJyrC1OUylwOnEOGYI7Px9PaRmJcIhk\nWxh3YQB3oJCc8nKrHG19PbHtR6y67l2WoXXWaY8erKO91qr9ngyFcebm4r/wAmJ19Rmz43v6uaU3\ncOnI6e49ei1TkzNNr6E+d+5cDMPgyJEjzJgxg/PPPx+n05l+/N///d97ffGKigpefvllAC666CJe\nfPHFY455+OGH0//+2c9+dkKNFxlsugZSvLmFeDAIxrGBE21osK4nNzeT7KyLflZROqC674fe9fXD\nu3eTbA1bFducDqtEqmkCJtHDhwlu3IR3aDkt2z6mbd8+AGJNTbgDAYZ/85s0b9lCPNhMPBQi1tiI\nmUxa1979fqvnDukd03A4wDBo/9vnpCJREuEQZiyOp7CIwv82haZNmzH/9jcitQfxDhua3sGtYEJl\nOoAToRBmAnKKAAAgAElEQVTJUIhUJIJvzBg8Zdakvs7Z8YFJk2ipqj5aXc80rW1VsTaDMUxrhj8c\nOyFOy9TkTNNrqN9555391Q6R09KpHr7t2nMM79ptVTrrCKSugdO10Ep4zx68xcV4SkuJ1NVlrEvv\n3harh/8hzjw/br+PRFs7DpcTR84QXHl+HA5HekZ8Itic8dxUWxvDb7qRvLHn0bRxE4ff/U8Ml9Pa\nDMXpsK6f5+fT/re/kUylMJwOEs0thHft6ij5aqRDP1JbS3DjZgyHYS1RKypMr13vPE9PSQlNGzbS\nXlNL7MhhPGedRXttLb6vfIW8sedlhHhGbztlpifuecvKMEzIHTmyxwlxWqYmZ5peQ/2SSy4BYNOm\nTRn3G4aB1+ulpaUlPflNxI5O9fBt1+B2+n0kQlbAQmbgdB1a959/PrFgU4/r0ru3pWsP3zt8OK7W\nUHpo2l1cTCoSwZmba12/LiggEWq1htE9HgomTUhfx3f4OqrBpV/ZxFNSgsPjsarNtbSQDIdIxOIk\nIxHcRYWkQiGMnCG48/PS26sC6fY4/b5jztM0wHC5cObm4vTlkjN8GO6C/IzzanjzreP+PA3DoHDK\n5ON+JlqmJmeaPl1Tf+yxx9i2bRuXXnoppmmyceNGhg8fTigU4u6771bJWLGtUz1827Xn6C0rw/OV\nQtwF+RmB0310oHja12mt3krDunUMGTYs3YPvqS3dS6l6v/IVnD4fsfp64qEQrtxca0vUw4dxFRaR\naGvDTKbwnDOGnDHnUL/6t9aQf00Njtxc/GPPI3DJFIIbN5EIh0lFoyRaWki2tlpD+h1135PJpLXR\nSjKJu7gYT1np0d3c6uvxNbcQDzZByipNm0omrSI4rSEMpwOH04WBgbesjJwu+8J3/5mBtVFLZx38\nLwpqzYKXM02fQt00TdasWZOeKFdfX8+SJUtYtWoVc+fOVaiLbZ3q4du+lDrtbXSga+U4b2lpxgzy\nRMeac3cgwJCRI/AUF1P/2hu0bN2KKxDAXRgg1dZOoqWF2OHDRA7UkGhrw+nxEKmtpfb5F/CedRaR\nmhrira04ojGcubkEN24i1tREsr2dRChkrUnvbHNHzXdCYQyfzyop2x7BW1iUruXetm8frZ/sINne\nzpCKCpo2bSa8czeR+jrCu3eTaG3F4fHgLgzgCRQeE9J54/6O0M5dtO3fR+6os8kfb+0SJyLH6lOo\nNzQ0ZMx8Lysro6GhAb/fb02WEbGpUz1825ee4/FGB7qHW15HaHadQZ4zfBg55eX4zzu3Y736TuKH\nDxOpPYjL7ye/cjzx5iCAdR0c0kPsiWAz+RddRHtNLalYHKffR7K9ndbaWkys9eUOlwt3IEA8GCQV\ni1nXzN0eHDk5uPPzcPr9uAsKcBfk07r1r9S/9jqtn31GpPYgBhA/fJh8xpNobsYEzHgCAwPD4SRn\n+HDcBfnW+vRuEwpjQWsiX7w5mN4lTkSO1adQnzBhAt/73veYMWMGqVSK119/nQkTJvCXv/yF3Nzc\nbLdRZMAMxPDt8UYHWrf+lXhzMCPc0jPIOwo3de4vHm1ooG3/PuBoaKcScUzTtIrJgLWELBrB4fUA\nkD/u7yicMhlMk0hDPYlQiEhNLYm2NlKxGInWVtyFRQypqMBdkE/7gRpceXnkVFRgOAxSHV8SXD5r\nT/RoQ4NV5rW1FbNjtnysKUjrto8p/Op/IxVuw5WfR8IAd14eyY7NYuDYCYVOvy+9Jl8z2EWOr0+h\n/rOf/YwXX3yRl156CafTyd///d/zzW9+k/feey9jSZqIfHnHGx3oqQfvKS7myMaNxINB6xr6sKGA\n9UUgd9TZBDdtxsAKds9ZZ2EADo8X3+jRGB4PidZWDKcT35gxnH3nHYS3fcyQUSNxFxYSazpC5ECt\nNaM9ZWImYhgGeIrPIu+iCygrLMKVn4e3rJRUIknDG38iGQ7jP/8r5I0fh1lVTTIStUI9kcBwOgAT\nw+3CyB0CDgeG04nnrLPwFBcTmDQpXQe+aeMmwnv24PL5cPp86S8rnecmIj3r8/D7FVdcwRVXXJFx\n3+WXX561homcqY43OtBTD97ExDDBHQjgHDKEnLIyCidPJr9yPHnjxxH86CNatv4VZ24upstF7MgR\nnDk5uIuK8JaXkztqZLoGe2cxmE6FU6YQraun9ePtYFrlWl1+v7U9qwmuPD+JFmtb1ERLKzlDh4IB\n8eZmWqu3dpwMODweSKYwHQ7chYX4vjKWeGMjpFL4zznnmL3Rm6uqiTTUpwvv5AwbRmDSpGMmFIrI\nsfoU6rfeemt6Mk88HufQoUNccMEF/Pa3v81q40TkqJ568I1v/TldpAWs9dqdXwgMw6D48qk4c4eQ\nDIVJRqOk2ttJRaPp8quFUyanX7/7SECssZHyGdeCw7Am1YXDeIcPI1pvFXyJ1tenr+MnQ9bWrd7y\nMqJ11h7v7vwCHF4PuWefbU2wi8VwF1mV5tp270mv0fdC+lp6+vw6N5wJhckpK0sHvoj0rk+h/vbb\nb2fc3rp1Ky+88EJWGiQiPeveg++6XakVqOXHDE13rQ8f3rUbT2EhOKytVL1lmeVXj7fNa+eytERL\nK678PNr3fw6GQXj3boD0hjKJUBjqjgZ9qj1CKhoDrEl2OcOGkXf++bgL8nEXBAh9+mn6enlg0qRj\n2tHZ7s791kXki53UNmjjxo1jyZIlp7otInICWqq3prcrTYTC+AOBY4amu/bu3QWBdFlayCy/6iku\nxjQMq1obEJgyOWOb14Iurxn88CPqX3u9o5RsK44cLyYmOWXlpBJxcoYPsyq9YeDp2McdrPXlBZVW\nedvgRx8R/vRTAIyOPd17arOG20VOTJ9C/f/+3/+bcXvXrl2cddZZWWmQiPRN1+1KIXMIu1NnKJsp\na7vU4MbNYFqhDXBk0yarUt2BGpxDhuC/4AJi9fUEN2221p87DGINjcesqTcNcAcKrA1f2iN4i87C\n5fenZ+Z3Kvpvl/Q4PyDW0Jhx2SDW0HhMm0XkxJ1UT33KlClcc801p7otInICTqQwTkv1VoKbt6Rr\npnfWYO8sPRtvChJvOrp+PREOZ2yW0rUITqyxMWPLVzi6oYorPw//eed+YS9bNdlFsqNPof6d73yH\nI0eOUF1dTTKZpLKykkAgkO22iUgvTmSYuqflcN7S0vRSMcPjtorDBIM4cnJw+HKJ1NSQikQB8JaX\np1+jayA7/T5r+LxDTllZn3rZGmIXyY4+hfp//ud/smTJEiorK0mlUvz4xz/mwQcf5Gtf+1q22ycy\noE71Lm2n0okMU3fvGXtKSjAxcebmEg8G8Z17LobDkR5yxzRJhMIYQHtHRTl3IEDDm2/hKSkhMHkS\nsYZGa4KbQcYQ/aluu4j0XZ9C/Ve/+hW/+c1vGDFiBAB/+9vf+M53vqNQF9s71bu09aQ/vjh07xmb\npklw8xZyhg61dl0sK6PwkinkjR+X3jzGf955YBgkw9YGME2bN+PK8+MtK6NoypT0+vbezqmzJr0r\nP4+cjtn2p8uXIhE76lOoJxKJdKADjBgxglQqlbVGiZwuTvUubT3pjy8O3XvG6e1MDet6eNf17d03\nj4kerCMebE7v0BY7fIRExzr344X08WrSZ+PcROSoPm11NGzYMJ577jlCoRChUIjnnnuO4V32Sxax\nq+4TuLIxoetUfnEwUymCH37EvqeeYd9TzxD88KMeN136ovPKG/d36SVwpmniHToU0zSJNjYS3rmL\nSE0tRzZtoqWqutdz6qkmvYhkT5966g8++CDLli3jySefxDRNvvrVr/Kzn/0s220TGXD9MaHrVM4E\nb6neSv1rr9NeWwtgzWDv4fr1F51X181jUu0RMAxcQ4YQiUSsjWAiEaL19ccN6c5zcvl8xJubcfp9\nX/rcROSL9SnU//3f/51/+Zd/yXZbRE47/TGhK79yvHWNe5M13G2aJqZpntS1586d0TolQ+Eeg/eL\nzqvrc7zl5WCaOIdYu7t1btmaDIWPG9KdXxIi9fX4OyrR5XSrYCcip16fQv2dd95h4cKFmuAikgWG\nYVi/Wx1ryINbtpz0lwlvaWm6dwzWkrOT6R1njB4Y1uYuAEc2biL0ySfEg0G8Y8eSN37ccc+peyU6\nEcm+PoV6IBDgqquu4qKLLsLr9abvX7FiRdYaJnImOVXX1fMrx2NiWpXjOFru9WRep7MdXYfnQzt3\n4fB6yT1nDE6/j9bqrZr4JnIa6VOoz5w5M9vtEDmjnarr6oZhEJgwgcCECV+qPccbKXAX5OM795z0\nbU18Ezm99DnU29raaG5u7nEmrYh8OYOlwprKu4qc3vq8ocuzzz5LYWEhhmGkJ/GsX78+2+0TOSMY\nhkH++HHpIjQtVdWnZaGWwfLlQ+RM1adQf/XVV3n77bcpLCzMdntEzlinqghNNivUqbyryOmtT6Fe\nWlpKXl7eCb94dXU1jzzyCKtWreKTTz5h2bJlOJ1OPB4PP//5zykuLs44fubMmfj9fgAqKio0EU/O\nKKdqslx/VKgTkdNTr6HeuY96fn4+s2fPZurUqTidzvTj3/nOd4773GeeeYY1a9YwZMgQwCpgc//9\n93PBBRfw4osv8swzz7B48eL08dFoFNM0WbVq1Zc6IZHB6lRdr47U1RE9WEciHMbl8xGpr9fSMpEz\nRK9lYl999VUAxo0bx9e+9rWMQP8iI0eOZOXKlenbv/zlL7ngggsASCaTGUvjAHbs2EF7ezvz589n\n3rx5VFVV9fm9ROwgv3I8hZMnkztyJIWTT24pGkCiNUR7bS3x5mbaa2tJtLSe4paKyOmq1556fn5+\nr73x3kyfPp0DBw6kb5d29Do+/PBDnn/+eV544YWM43NycliwYAE33XQT+/bt49vf/jZr167F5erT\nFQKRQe9UXa925eeRM3wYyVAYp9+HK//EL52JyODUa2Ke6pm3b7zxBk888QRPP/00RUVFGY+NHj2a\nUaNGYRgGo0ePJhAI0NjYyNChQ09pG0TsLqesLL0jWudtETkz9BrqO3fu5Morrzzm/pNZ0vaHP/yB\nl156iVWrVhEIBI55fPXq1Xz22WcsXbqU+vp6QqEQJSUlfX59EbGc7LKz/tjXXUSyq9dQHzVqFE8/\n/fSXfpNkMsmDDz7I0KFDufPOOwGYMmUKd911F4sWLWLhwoXMmjWLxYsXM2fOHAzDYPny5Rp6FzkJ\nJzuMr1nzIoNfr6npdru/1L7pFRUVvPzyywBs3Lixx2Mefvjh9L8fffTRk34vEflyTuW+7iIyMHqd\n/T5x4sT+aoeIDLDuS+hUAlZk8Om1p/7jH/+4v9ohIgNMJWBFBj9dtBYRQCVgReyg1+F3ERERGTwU\n6iIiIjahUBcREbEJhbqIiIhNKNRFRERsQqEuIiJiEwp1ERERm1Coi4iI2IRCXURExCYU6iIiIjah\nUBcREbEJhbqIiIhNKNRFRERsQqEuIiJiEwp1ERERm1Coi4iI2IRCXURExCYU6iIiIjbhGugGiIjI\n6cVMpWip3kq0oQFvaSn5leMxDGOgmyV9oFAXEZEMLdVbadq8GYC2zz8HoGBC5UA2SfpIw+8iIpIh\n2tDQ6205fSnURUQkg7e0tNfbcvrKaqhXV1czd+5cAPbv38+cOXP41re+xU9+8hNSqVTGsalUih//\n+MfMnj2buXPnsn///mw2TUREjiO/cjyFkyeTO3IkhZMnk185fqCbJH2UtVB/5pln+NGPfkQ0GgVg\nxYoVLFy4kN/85jeYpsn69eszjl+3bh2xWIyXXnqJ733vezz00EPZapqIiPTCMAwKJlRSOv0bFEyo\n1CS5QSRroT5y5EhWrlyZvv3xxx9zySWXADB16lTef//9jOO3bNnCZZddBkBlZSXbtm3LVtNERERs\nKWuhPn36dFyuo5PrTdNMf9vz+Xy0trZmHB8KhfD7/enbTqeTRCKRreaJiIjYTr9NlHM4jr5VOBwm\nPz8/43G/3084HE7fTqVSGV8KREREpHf9FuoXXnghGzZsAODdd99l8uTJGY9PnDiRd999F4CqqirG\njh3bX00TERGxhX4L9R/84AesXLmS2bNnE4/HmT59OgCLFi2itraWadOm4fF4uPnmm1mxYgWLFy/u\nr6aJiIjYgmGapjnQjThZBw4c4Morr2T9+vVUVFQMdHNERESy6otyT8VnREREbEKhLiIiYhOaXi5i\nY9ptS+TMolAXsTHttiVyZtHwu4iNabctkTOLQl3ExrTblsiZRcPvIjbWubtW12vqImJfCnURG+vc\nbUtEzgwafhcREbEJhbqIiIhNKNRFRERsQqEuIiJiEwp1ERERm1Coi4iI2IRCXURExCYU6iIiIjah\nUBcREbEJhbqIiIhNKNRFRERsQqEuIiJiEwp1ERERm1Coi4iI2IRCXURExCYU6iIiIjahUBcREbEJ\nhbqIiIhNuPrzzV599VV+97vfARCNRvnkk0947733yM/PB+C5557jlVdeoaioCICf/vSnjBkzpj+b\nKCIiMmj1a6jfcMMN3HDDDYAV2DfeeGM60AG2bdvGz3/+cy6++OL+bJaIiIgtDMjw+1//+ld27drF\n7NmzM+7/+OOPefrpp5kzZw5PPfXUQDRNRERk0BqQUH/qqae44447jrn/mmuuYenSpfz6179my5Yt\nvPPOOwPQOhERkcGp30O9paWFvXv38tWvfjXjftM0+cd//EeKiorweDxcfvnlbN++vb+bJyIiMmj1\ne6hv2rSJSy+99Jj7Q6EQ1157LeFwGNM02bBhg66ti4iInIB+nSgHsHfvXioqKtK3//jHP9LW1sbs\n2bO55557mDdvHh6Ph0svvZTLL7+8v5snIiIyaPV7qP/zP/9zxu0ZM2ak/3399ddz/fXX93eTRERE\nbEHFZ0RERGxCoS4iImITCnURERGbUKiLiIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVE\nRGxCoS4iImITCnURERGbUKiLiIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGxCoS4i\nImITCnURERGbUKiLiIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0o1EVERGzC1d9vOHPmTPx+\nPwAVFRWsWLEi/djbb7/NY489hsvl4sYbb+Sb3/xmfzdPRERk0OrXUI9Go5imyapVq455LB6Ps2LF\nClavXs2QIUOYM2cOV1xxBcXFxf3ZRBERkUGrX4ffd+zYQXt7O/Pnz2fevHlUVVWlH9u9ezcjR46k\noKAAj8fDpEmT2LRpU382T0REZFDr1556Tk4OCxYs4KabbmLfvn18+9vfZu3atbhcLkKhEHl5eelj\nfT4foVCoP5snIiIyqPVrqI8ePZpRo0ZhGAajR48mEAjQ2NjI0KFD8fv9hMPh9LHhcDgj5EVERKR3\n/Tr8vnr1ah566CEA6uvrCYVClJSUAHDOOeewf/9+gsEgsViMzZs3M2HChP5snoiIyKDWrz31WbNm\nsXjxYubMmYNhGCxfvpw//elPtLW1MXv2bO677z4WLFiAaZrceOONlJWV9WfzREREBrV+DXWPx8Oj\njz6acd/EiRPT/77iiiu44oor+rNJIiIitqHiMyIiIjahUBcREbEJhbqIiIhNKNRFRERsQqEuIiJi\nE/2+oYuIiEi2mKkULdVbiTY04C0tJb9yPIZhDHSz+o1CXUREbKOleitNmzcD0Pb55wAUTKgcyCb1\nKw2/i4iIbUQbGnq9bXcKdRERsQ1vaWmvt+1Ow+8iImIb+ZXjATKuqZ9JFOoiImIbhmGcUdfQu9Pw\nu4iIiE0o1EVERGxCoS4iImITCnURERGbUKiLiIjYhEJdRETEJhTqIiIiNqFQFxERsQmFuoiIiE0M\n6opyyWQSgLq6ugFuiYiISPZ15l1n/nU3qEO9sbERgFtuuWWAWyIiItJ/GhsbGTVq1DH3G6ZpmgPQ\nnlMiEomwbds2SkpKcDqdA90cERGRrEomkzQ2NnLxxReTk5NzzOODOtRFRETkKE2UExERsQmFuoiI\niE0o1EVERGxCoS4iImITg3pJWzbMnDkTv98PQEVFBStWrEg/9vbbb/PYY4/hcrm48cYb+eY3vzlQ\nzTxhr776Kr/73e8AiEajfPLJJ7z33nvk5+cD8Nxzz/HKK69QVFQEwE9/+lPGjBkzYO3ti+rqah55\n5BFWrVrF/v37ue+++zAMg/POO4+f/OQnOBxHv7OmUimWLl3Kp59+isfj4YEHHuhxOcjpouu5ffLJ\nJyxbtgyn04nH4+HnP/85xcXFGcf39v/2dNL1vLZv385tt93G2WefDcCcOXO4+uqr08cO5s/snnvu\n4dChQwDU1NQwfvx4fvWrX2Ucf7p/ZvF4nCVLllBTU0MsFuP222/n3HPPtcXvWU/nNmzYMHv8npmS\nFolEzOuuu67Hx2KxmPn1r3/dDAaDZjQaNW+44QazsbGxn1t4aixdutR88cUXM+773ve+Z/71r38d\noBaduKefftq89tprzZtuusk0TdO87bbbzP/6r/8yTdM077//fvOtt97KOP7NN980f/CDH5imaZof\nffSR+b/+1//q3wafgO7ndsstt5jbt283TdM0/+M//sNcvnx5xvG9/b89nXQ/r5dfftl89tlnj3v8\nYP7MOgWDQfMf/uEfzPr6+oz7B8Nntnr1avOBBx4wTdM0m5qazMsvv9w2v2c9nZtdfs80/N7Fjh07\naG9vZ/78+cybN4+qqqr0Y7t372bkyJEUFBTg8XiYNGkSmzZtGsDWnpy//vWv7Nq1i9mzZ2fc//HH\nH/P0008zZ84cnnrqqQFqXd+NHDmSlStXpm9//PHHXHLJJQBMnTqV999/P+P4LVu2cNlllwFQWVnJ\ntm3b+q+xJ6j7uf3yl7/kggsuAKw1ql6vN+P43v7fnk66n9e2bdv4y1/+wi233MKSJUsIhUIZxw/m\nz6zTypUrufXWWyktLc24fzB8ZldddRV33303AKZp4nQ6bfN71tO52eX3TKHeRU5ODgsWLODZZ5/l\npz/9Kd///vdJJBIAhEIh8vLy0sf6fL5j/ggNBk899RR33HHHMfdfc801LF26lF//+tds2bKFd955\nZwBa13fTp0/H5Tp69cg0TQzDAKzPprW1NeP4UCiUHjYDcDqd6c/2dNP93DoD4cMPP+T555/nf/7P\n/5lxfG//b08n3c9r3LhxLFq0iBdeeIERI0bw2GOPZRw/mD8zgMOHD/PBBx9www03HHP8YPjMfD4f\nfr+fUCjEXXfdxcKFC23ze9bTudnl90yh3sXo0aP5h3/4BwzDYPTo0QQCgXQpWr/fTzgcTh8bDocz\nQn4waGlpYe/evXz1q1/NuN80Tf7xH/+RoqIiPB4Pl19+Odu3bx+gVp6crtf1wuFweq5Ap+6fXyqV\nOuaP8OnsjTfe4Cc/+QlPP/10et5Dp97+357Opk2bxsUXX5z+d/f/c4P9M1u7di3XXnttj9UuB8tn\ndvDgQebNm8d1113HjBkzbPV71v3cwB6/Zwr1LlavXs1DDz0EQH19PaFQiJKSEgDOOecc9u/fTzAY\nJBaLsXnzZiZMmDCQzT1hmzZt4tJLLz3m/lAoxLXXXks4HMY0TTZs2JD+YztYXHjhhWzYsAGAd999\nl8mTJ2c8PnHiRN59910AqqqqGDt2bL+38WT94Q9/4Pnnn2fVqlWMGDHimMd7+397OluwYAFbt24F\n4IMPPuCiiy7KeHwwf2ZgndPUqVN7fGwwfGaHDh1i/vz53HvvvcyaNQuwz+9ZT+dml9+z0/Mr1ACZ\nNWsWixcvZs6cORiGwfLly/nTn/5EW1sbs2fP5r777mPBggWYpsmNN95IWVnZQDf5hOzdu5eKior0\n7T/+8Y/pc7vnnnuYN28eHo+HSy+9lMsvv3wAW3rifvCDH3D//ffzy1/+kjFjxjB9+nQAFi1axMKF\nC5k2bRrvvfceN998M6Zpsnz58gFucd8kk0kefPBBhg4dyp133gnAlClTuOuuu9Ln1tP/29O1d9TV\n0qVLWbZsGW63m+LiYpYtWwYM/s+s0969e48Jh8H0mT355JO0tLTw+OOP8/jjjwPwwx/+kAceeGDQ\n/551P7dkMsnOnTsZNmzYoP89U+13ERERm9Dwu4iIiE0o1EVERGxCoS4iImITCnURERGbUKiLiIjY\nhEJd5AyzYcMG5s6d2+fjT+RYERlYCnUR6dXGjRsHugki0kcKdREhkUjwox/9iNmzZ3PllVfyz//8\nz0QiER544AEAbrrpJsCqIjZr1iyuv/56vvOd79DU1ATAFVdcwb/8y78wa9YsrrnmmvRGHp988gk3\n3XQTM2bM4NZbb6Wuro57772Xl156Kf3ec+fOpbq6up/PWMSeFOoiwkcffYTb7eall17iz3/+M9Fo\nlP/3//4fP/rRjwB45ZVXOHLkCI8++ijPPvssv//97/nv//2/88gjj6RfIxAIsHr1am6++eb0Tn/f\n//73+d//+3/zxz/+kauvvppf//rX3HjjjaxZswaw9ho/cuQI48eP7/+TFrGh06/GnYj0uylTphAI\nBHjhhRfYs2cP+/bto62tLeOY6urq9CYYYG3WUVBQkH68c8vN8847j7feeosjR47Q2NjI1772NQC+\n9a1vAdYGQvfffz8HDhzgD3/4A9ddd11/nKLIGUGhLiKsX7+e//N//g/z5s3jhhtuoKmpie4VpJPJ\nJBMnTuTJJ58EIBqNZuzI1bn/dOfWnG63O+P50WiUhoYGRowYwfXXX8/rr7/O2rVr+bd/+7dsnprI\nGUXD7yLCBx98wP/4H/+DG2+8keLiYjZt2kQymQSO7ok9fvx4qqqq2Lt3LwCPP/44Dz/88HFfMy8v\nj/Lyct577z3A2gXrX//1XwG44YYbePHFFykvLx90GyOJnM7UUxc5A3XfOnjcuHFs2LCBtWvX4vF4\nqKys5MCBAwBceeWVXHfddbz66qssX76chQsXkkqlKCsr4xe/+EWv7/OLX/yCpUuX8vDDD1NYWJj+\nEgWJCicAAAB1SURBVDB06FDKy8uZOXNm9k5S5AykXdpEpF+ZpklDQwNz587ltddew+PxDHSTRGxD\nw+8i0q/efPNNrrvuOr773e8q0EVOMfXURUREbEI9dREREZtQqIuIiNiEQl1ERMQmFOoiIiI2oVAX\nERGxCYW6iIiITfz/AAU0nn7Zw0sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb9f0b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.regplot('Latency', 'Throughput',\n",
    "           data=pd.DataFrame(X, columns=['Latency', 'Throughput']), \n",
    "           fit_reg=False,\n",
    "           scatter_kws={\"s\":20,\n",
    "                        \"alpha\":0.5})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# estimate multivariate Gaussian parameters $\\mu$ and $\\sigma^2$\n",
    "> according to data, X1, and X2 is not independent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 14.11222578  14.99771051] \n",
      "\n",
      "[[ 1.83862041 -0.22786456]\n",
      " [-0.22786456  1.71533273]]\n"
     ]
    }
   ],
   "source": [
    "mu = X.mean(axis=0)\n",
    "print(mu, '\\n')\n",
    "\n",
    "cov = np.cov(X.T)\n",
    "print(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0, 0],\n",
       "        [0, 1],\n",
       "        [0, 2]],\n",
       "\n",
       "       [[1, 0],\n",
       "        [1, 1],\n",
       "        [1, 2]],\n",
       "\n",
       "       [[2, 0],\n",
       "        [2, 1],\n",
       "        [2, 2]]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example of creating 2d grid to calculate probability density\n",
    "np.dstack(np.mgrid[0:3,0:3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFbCAYAAAAuggz2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0W/Wd9/H3lWTJsmRFdrzHSeyYhGyYENKwlKYtNA1D\nh0mgaRlKw3MapoeHKS05LbSQhqUNhbbQ6cLA0GF6zkxpeqBD+9DMzCm0DaSZEggESEI2yOKQxYu8\n25JtSVe6zx+OhZ04iRIs29f+vP5JLFu6X8lX/uj3u7/FsCzLQkRERGzJMdIFiIiIyLlTkIuIiNiY\nglxERMTGFOQiIiI2piAXERGxMQW5iIiIjbky9cCJRII1a9ZQU1ODYRh85zvfwePxcPfdd2MYBtOn\nT+f+++/H4dBnCRERkXOVsSB/+eWXAXjmmWfYsmULP/7xj7Esi1WrVnHJJZdw3333sWHDBhYvXpyp\nEkRERMa8jDWHP/WpT7F27VoAamtrCQQC7Nq1i4ULFwKwaNEiNm/enKnDi4iIjAsZa5EDuFwuvvWt\nb/GnP/2Jn/3sZ7zyyisYhgGAz+ejs7PzlPft6elh586dFBYW4nQ6M1mmiIjIqJBIJGhsbGTu3Llk\nZ2endZ+MBjnAD37wA+68804+//nPE41GU7dHIhECgcAp77dz505uuummTJcnIiIy6qxbt44FCxak\n9bMZC/Lnn3+ehoYGbr31VrxeL4ZhMHfuXLZs2cIll1zCpk2buPTSS095/8LCQgD++an/oKi4JFNl\nioiIjBqhhnpu//L/SWVgOjIW5J/+9Ke55557uOmmmzBNk9WrV1NVVcW9997LP/3TPzFt2jSWLFly\nyvv3dacXFZdQWjYpU2WKiIiMOmdzSTljQZ6Tk8NPf/rTk27/1a9+lalDioiIjDuaxC0iImJjCnIR\nEREbU5CLiIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CL\niIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CLiIjYmIJc\nRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CLiIjYmIJcRETExhTk\nIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CLiIjYmIJcRETExhTkIiIiNqYg\nFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbc2XqgePxOKtXr+bYsWPEYjFuu+02SktLufXW\nW6moqADgxhtv5JprrslUCSIiImNexoJ8/fr1BINBHnnkEdra2li2bBlf+cpX+NKXvsTKlSszdVgR\nEZFxJWNBfvXVV7NkyRIALMvC6XSyc+dOampq2LBhA1OnTmX16tX4/f5MlSAiIjLmZewauc/nw+/3\nEw6H+drXvsaqVauorq7mm9/8JuvWrWPy5Mk8/vjjmTq8iIjIuJDRwW51dXXcfPPNLF26lGuvvZbF\nixczd+5cABYvXszu3bszeXgREZExL2NB3tTUxMqVK7nrrrtYvnw5ALfccgs7duwA4NVXX2XOnDmZ\nOryIiMi4kLFr5E8++SQdHR088cQTPPHEEwDcfffdPPTQQ2RlZVFQUMDatWszdXgREZFxIWNBvmbN\nGtasWXPS7c8880ymDikiIjLuaEEYm7JMc6RLEBGRUSBjLXLJnMibW4nX1ZJVWobv4gUjXY6IiIwg\ntchtxjJN4nW1AMTratUyFxEZ5xTkNmO4XGSVlgGQVVqG4VKniojIeKYUsCHfxQuwTFMhLiIiapHb\nlUJcRERAQS4iImJrCnIREREbU5CLiIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxM\nQS4iImJjCnIREREbU5CLiIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJj\nCnIREREbU5CLiIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREb\nU5CLiIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CLiIjY\nmIJcRETExlyZeuB4PM7q1as5duwYsViM2267jfPOO4+7774bwzCYPn06999/Pw6HPkuIiIicq4wF\n+fr16wkGgzzyyCO0tbWxbNkyZs6cyapVq7jkkku477772LBhA4sXL85UCSIiImNexprDV199NXfc\ncQcAlmXhdDrZtWsXCxcuBGDRokVs3rw5U4cXEREZFzIW5D6fD7/fTzgc5mtf+xqrVq3CsiwMw0h9\nv7OzM1OHFxERGRcyeoG6rq6Om2++maVLl3LttdcOuB4eiUQIBAKZPLyIiMiYl7Egb2pqYuXKldx1\n110sX74cgNmzZ7NlyxYANm3axIIFCzJ1eBERkXEhY0H+5JNP0tHRwRNPPMGKFStYsWIFq1at4rHH\nHuOGG24gHo+zZMmSTB1eRERkXMjYqPU1a9awZs2ak27/1a9+lalDioiIjDuaxC0iImJjCnIREREb\nU5CLiIjYmIJcRETExhTkIiIiNqYgFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CLiIjY\nmIJcRETExhTkIiIiNqYgFxERsTEFuYiMSZZpjnQJIsMiY/uRi8joYpkmhmt8vOUjb24lXldLVmkZ\nvosXjHQ5Ihk1Pt7VIuPceAo2yzSJ19UCEK+rHVcfYGR8Ute6yBg3WLCNZYbLRVZpGQBZpWUKcRnz\ndIaLjHF9wdbXIh8Pwea7eIFa4jJupHWW79u3j+nTpw+4bdu2bcybNy8jRYnI0BqPwTaenquMb6c9\n0998802SySRr1qzhe9/7HpZlAWCaJg888AAvvvjisBQpIh842tL9Ie4dpzzfO2S1iMjIO22Qb968\nmddff51QKMRPf/rTD+7kcnHDDTdkvDiR8ezDBfYHTmyJD/a4CncR+zptkH/1q18F4Pnnn2fZsmXD\nUpDIeDZU4d3H3PEWVqgeo6gEV/X8tI/bP9jHW5e8iN2k9e7csmULW7ZsOen2hx9+eMgLEhlvhjq8\n+1imiRWq7/1/qP6sArmvJnPHWwTDzeNi2pqIXaX1rl64cGHq/6ZpsmHDBqZNm5axokTGg7MN8LNt\nGRsuF0ZRSapFfrat6r4PAq0AB97HNfl8JhflntVjiEjmpfXOvu666wZ8vXz5cm688caMFCQy1p1L\nC/xMXeR7mjoGv2NxBRRXMKs4/6yPOdgHgb7ay/O96nIXGSXO6V144MABQqHQUNciMuadLsRPFYx9\nLeO6cA+EDxEprgDnmd+6nr07cIXqMYtK2OOs/uAbCXPA/WcVBE75GK7q+SR7enBkZw+4/dDGV7BC\n9RRVTVWXu8gISyvIZ86ciWEYqeln+fn5fP3rX89oYSJjyZla4adqcfe1tD3ZAVzhHsyCorRCnIRJ\n9s63cIY7cYVqiU6fDU5Xb7g3hTALiojOrB5wDDg51PvqSh6vq29VuL5r7yF1uYuMuLSCfO/evZmu\nQ2TMOlOIDzYobW9b14Cfic6sJnpCS3pQp/qZRG8Au5pCqX8He7xUqCdMZk4MDKjLfHsrVnOot5u9\nr8t9YlHqOWoKm8jISCvIOzo6eOyxx3jttddwuVwsWrSI2267jewTuttEZKB0rof3vxZ9JDtA9IQQ\nTzlNiL9a00rRoT34OluIBAtpqJhNcdEs/K56smJRYi/8makzpmAWFKVa5KnHOx7yfV/3tdprCoog\nO8Dkng6YWIjV3PshwArV4/rEp4nHYiS3b8XcvhXX/IUcPd6ToEAXGV5pBfldd93FtGnTePTRR7Es\ni9/+9rd8+9vf5kc/+lGm6xOxrbMZ1Lav7LzegWlpdJu/WtOK1e/nig/tZt6+bQSaamkvLAfASJg0\nVMymsayKit2vAfD+e4epueCjMGly7/0Pd/LJzv1k73obHA565s4nOn02ruOtcFdTiMhln2A/gNPF\ndKcz1f2feGcbic0bsdracOTlk6yrxZpdnXrefWGuAXH2od+VfaX1Wzt27Bg///nPU19/+9vf5m//\n9m8zVpSI3Z1NiKe6s8/Q4racLooP7aayrTHV6jYSJv6Wetw9EcDC09VOJHdOKuiTbg+RYCG+4/ex\nnC6M4y3wkoPvkHztD1jhdrpz8zjqymdKIoFn7w4c3V30zKoeUFPfh43zJ+SQ+OtLGIEgtLWR9Obg\nKi0jsXtHKuiPzK4m/8i7xI4dxT2pXAPiRrnxtM3tWJRWkE+dOpWtW7eyYEHvL3jv3r1MnTo1o4WJ\njAennDZ23ObDnanw7grkk9PRAoCvrRFHLErS7SGcX0JuS4io10/9tLmEps4CelvlltOVCnzobb37\n2hqJ5Obja28Cw+g9kGXRFcjj8L7DFMQdFPoDWNk5J11z9+x9h6OtzRhdYUrKp5AMBDF8fqxEApob\nATDfep1E7WFqd+7A8GSTe/QIORfOU2tvlNL+7faX1m/r8OHDfPGLX6SyshKn00lNTQ0TJkzgyiuv\nxDAMNmzYkOk6RWwj3dZ4OiFuJEx8bb0BmdPR0hvAnS24e7qo2P0akWAhf7ZKcFQWAlDZcIC8g38k\nO9ZNj9tLqy+fmqIqKkMHyO9sIquzEaqm4W/vbaG3F5bTUTCJ2mkX0FB1AcWHdlN8aA+d4W4actvI\n7x/iO7bie+tVEv5czJJy9p8/j/Pe3UYyHsdoboTgRKzmRgzAMFxYXV3gdhPuMSk6t5dShsF43OZ2\nrEnrN/bkk09mug6RMWEoQnzz4c7U/y2na0DXeEPFbByxKFmbNkA4Ak0RHJWFJB1OXGaMvEgLRjJJ\ncXsdhwsqmFG3h4KOEBPDTYQmlIAB8QMHsAyDg6U5bKj8FHMr8lJd8Y1lVfhbGnqP7XDwak0rAJ/o\nOoTvrVchFsUZhp68ieD107RzJ4GGoxi5AbIWXo5RVIyjtAwrVI+zej5GjhdHySSOdcQpz1dAjFbj\ncZvbsSSt31pRURGbN2+mtbV1wO3aSEXk7J0uxF+taT3pWnn/rvF3anq71qt8+eRFWmj15ZN0OKmq\nf4+8SEuqJR4KFIMFRjJJfriZovYGHEmTvWVzyYu0YDkc5EVaqKGK3QebMF1uqurfY643Slash1h2\nDpFgIUVH3sPf0kBrSwPZgVycHRaR+ZcRnXMRxHpwNYfoiZtkh+pJmjEczY04P/FpmF2N4XIpHGxE\nvyf7Sus3d8cdd9DY2EhVVRVG3zU1FOQiQ+nApi0DBrL1F3l9K3mRFqp8+RwomcGBkhk4kgmSDifu\nWDd5kd6A73F76cjuXdRlUtP75HW14o11E3O6KWs9Rkd2kGMFU1IfAj66dyNFHQ2Ecovp8Xg5FnEw\nqcDHodmXYjmdVL7zCpbDAVjUtUcp9fs+KMrpIunNwRHtocvrY4JlpJZy7Vs45sTtUzU1TWTopRXk\nBw8e5IUXXsh0LSK2lk63ev8FVwa0vI9fCzeSiQED2QB2HWhk/vGgzou04DJjmC43AB/bvYGijgZa\nOk0O+kpozApQHD1EGxZzWo9SnzWBCXEDr9WNRZLKmu281e2l5uJLmV67h9nHdhJ1ecjt7qQxUEjS\ncLIr53wOHulgTlVhqlu/ftpcSva/g7lnD8G97xCpPUz46uuJLPo0rvpajFgP3d1xymDAwjEnrguv\nMBcZemkF+ZQpU6itraWsrCzT9YiMeSctk5ow2Xysm5k9XUxoPErCmQWQGsiGw0nr8a707Fg3Fx3a\niiceJRhppvzoPlqzfLiz/GzOm8GnGt/hsra9YEHCcFAcb+dwdiEhd4AFHQdozsqlINZJ09YtXNj4\nJt5kFx5HlFCgiIZgGUXtdcyo28v0unfZ13A+DZdfjrO7i6TbTcn+bWRFI3Q7Xbhrj0J3uHfFuWkz\n8L3+VwBCf91EMNaFI5iPk8HnJivMRYbWaYN8xYoVGIZBS0sL1157LTNnzsTpdKa+/8tf/jLjBYqM\nFXuaOiBhDlgmlV1v42ptpjieQyw7h+bSSibW1WAkk70t89zegWwHSmbgjnVz1c4/UtIWprC9ga5Y\nEm8iSlbSpMnlB2BKd4iE0dsVHnJP4F1fCQmnm40T51KbncfMSC2Tepop62nBb3aTSJgczp7ArqpL\nyQ83gQXeeBe53Z34omHiR7ZTmOOivbCcuhkXk9dwBE93hPbmVnxvvpb6MGIWFOEK1eFsayHS3UVO\nezuuCxfouqvIMDjtu+yrX/3qcNUhMj44XR8sk5o3EVdrMwdaevDRk5on3l5YjuVwsLPbQ3JC7wfn\nqvr3yA83MaGrjYTDQTSaIOp00+30EDaymN5dz/99/0XCTi/dDjfdjiwO5RSTcLppcAf527otTOlp\nZmK8nUZPHgXRdnLN3ksBSYeL2mNtmLEwnrwAYJDb3Um328uklqO0mvlM4Cj7Lvw4h+dcRqDhMPmh\nIzS+u59CIJowUy1zV6gOx7HDGFinfRnUKhcZOqcN8oULFwLwxhtvDLjdMAw8Hg8dHR0EAqfeAlFE\nTjZgA5S9O6DlMJFgIaHJMwBSq68dONx7Pd2RTPSONDcc+HrCmJFu2rJ8vB2owBfv5uLOGrAM8uNh\n9vjLeTswlaDZTb0nnz255Sytf43FzdsxkkmcJEkYLjyJGFGHCwyDFlcORbE2ynpayDnaw77Zl/Je\n6UyC3W3kh1vIi7TQ4Col4c0hK9pFfugISYcTT08X75g5TOu71u/OpmfWhTg72mkLTKCwOTRo17pG\nsosMrbTeTY8//jg7d+7ksssuw7IsXn/9dSZNmkQ4HOaOO+445XKt27dv59FHH+Xpp59m9+7d3Hrr\nrVRUVABw4403cs011wzZExEZbh8qkI6HX3T6bGo8kyk68h6V77xCVyCf+mkXYDldqVHpAK2+fCZ2\nhkg4XdRn5xE0u3BicTC3jEIzzPmRY7gTcWZ1HsWbjHHEW4hlOHC3R7movQbDsvAmY3QbWeQkurEM\nBzHDTdjp4e1gFRVdIaZ31dGc5WdCtIO3Kj+CI5nASCZxWCYJpxtHLErM66e5bBqengjvz/wI9dMu\nYFq/pxWdeQEArtbm1Aj2/vpv18onPnpur52IDJDWXyHLsli/fn1qsFtDQwOrV6/m6aefZsWKFYMG\n+VNPPcX69evxenu7z3bt2sWXvvQlVq5cOYTli4yMoVibum/QW1E8B19nC4GmWgqPvIdlAU4nnppD\nA1Zo21p1GW4zjv/9A3Q5PRTGOqmK1NPp8rLTP4WqSD1eswfDMAjGOtjtn0Se2U2zO5dgPEyW4cIw\nINeMcjCnhJjDyV/zZrPXX05RtJ19OWXkJrrZ3QLJKieVoQNUNewjO97N/pLz2XosQmGwdwW5SG4+\nDVUXDPp8zIIiIpd9Aldx/oDvn7hd65FQp/YxFxkCjnR+KBQKDRixXlxcTCgUwu/3Y1mDXwubMmUK\njz32WOrrnTt3snHjRm666SZWr15NOBz+kKWLjIzB1qY+a/0Gvfk6W+j2T8DTHSaa7cPf3oi/pSG1\nQpthJVPzxOvyJvGuv4wD3mICZhcdWTlgQCDeRdzppDTeTllPC9O6Qlx6+HUc9Ud4vTubLck8ui3I\nScbwmd0URtvodGYzuaeZG2pfYW7nYSwD/jd/FjsmVOJIJijoCBF3ZRHOziXmzMKRTNBQMZuaCz46\nIMQ3H+48eRAfJy9807ddKzBoa11Ezk1aQX7RRRfxjW98g40bN/LSSy/xjW98g4suuoiNGzeSk5Mz\n6H2WLFmCq98btbq6mm9+85usW7eOyZMn8/jjjw/NMxAZZn1rUwPnvja104WZNxHonWZWd9482gvK\n8PREyIpFCecXYzkcqRXa2j29Y1HyIi0kHS4iWV5eyTuf2uyJlPS0MaWrgbx4hIasCZiGAwdJpiQ7\nqUq2M8NsZbIVwYuJBTix8JtdgMHUrgYC8Qhxh4t69wR25fZuhlQZOkB+pJksM47LjDMx3Exl6ADA\ngC1UoXdzlr5BfMDAvc5P4Kqej+sTnz5pfrmInLu0/gJ997vf5ZlnnuHZZ5/F6XRy+eWX8/nPf55X\nXnmFH/7wh2kdaPHixamBcYsXL2bt2rXnXrXICDuXtalnFQRSrVTP8WlnZt5EquZcRKimlZjXT9Ok\n83oDfPIMuo61UdjRwPm1e0g4nBR3NlCXNwnooMEdZE9uOd5EN8vqX8NjxUmacCy3Apdl4rYSeIkS\nNKPkEWO/YwITE934jCSGw0Gn2095TzNdTg8WBs1ZudRlF5AwnLiTMfI7m2gMFGMkEvhiXWD0foiI\nHd9RrU/fbmqe6JSBg/hOXPCmH7XERYZWWu+oUCjElVdeyZVXXjngto9//ONpH+iWW27h3nvvpbq6\nmldffZU5c+acfbUio8iJgVSe7x10dbcTA9+z6218b/yVhD8Xkgmix8Ox/+YoAPnhJnxdHRR2NIBh\n4I118dr0j3KwzcuMcB1XNm6j0TOBuOEgx0qSNJy4zThN7gn4zW5yzW4mueJ0Ody4vbm868jHBILH\nr5sbWBz1FuK0EmycOJeYw83cjve5qOMgU91ROnICWA4n4exccqIR2rxBvP3C2RGLpnZmczWFUiHe\n/1o5V1yRgVdeRPpLK8i/+MUvptZYj8fjNDU1MWvWLH7729+mfaAHHniAtWvXkpWVRUFBgVrkYhvn\nOjrdMk0Su3ekRmm7quf3XktubSbhzyXr2GFIJvHs2w05lanNUSyni+JDu6mKHMHRGcKwLCwMvLFu\nLn/3L1zk9jG9/nWyrATNWX7+kj+Xj7TvoyUrl+k9dfQYWUyMdxJxevAlonS6A2wNTmNGpJ6A2U1L\nlo93AlOZ23mY6o5DNLlzmRGuY09uOSXRFgLxCKY3m6bcQlr8BQS72mjLCXKweDp9H7/7WuLuni5i\n2TkfdKefcK1cU81EMi+td9hLL7004OsdO3awbt26M96vvLyc3/zmNwDMmTOHZ5555hxKFBk55zo6\n3dzxFsn6YyRrj+Esn4IVqscyTWYV53OwoAiSCUhamGWTcYXqufyjs9l8rDs1h9zfUk8iy01i0hQO\ntQTx93QS9gYo7GzkaL6HAiNGe9JJTiLGH4ovpjY7jwlmN7ktPeQmukkaDhxWEsOymBjr5FONO+hy\nZuNNxggSwepw9t4n3kWdJ4/iWBt7KMebiDMp2oKnw0V9sIz9pTNxJBNUhg4wv+YNgs4KQpNnpFri\nsewcDs2+lJLzCnqfeP8FbwqKFOIiwyCtwW4nqq6uZteuXUNdi8ioci6j08vzvalpVobDiQFYieSA\nUdrR6bOJfPQqeqovxlV/FFeoFs++3Vw+5YOpWOH8EqLeXKI5Abac/zH2ls8h6s4mYTjJj7QScefg\n9Hl4J7eCJaFtLGw/gM/s5v+VXsIhbyGvTziPlqxcwGJivJPJ3U04rCSdrmwOewuwjm9i+L63d4Ba\ngzvYW5vTRUdBKaHcQqIuD45kAiA1ar4vwPu6/yPBQi7tC/HjojOriVz2id515EUk49L6uPzP//zP\nA77ev38/EydOzEhBIqNF3+j0vhZ5uq3LvmlWVqge5/yFOI/vzQ29LfXzQvUcyQ4QnT4bV6gO+OAa\nc1+XdSQ3n+2fXA6AP5GgZ2OIbncO3miEvHALudEw7ngPM0sm4jAsnJ0RSqKtvO8t5qg7j/JYG4FE\nN9lWAgfgtkwmxjtoMiYA8Hagij255VzYfoD8eASAWZ1HWWiFKGmtI5rloSFv0oAFaeZ6o3QF8rGc\nrgGXAaoGexFOMdBNRIbeOb3bPvKRj/CZz3xmqGsRGXXOZXR6eb6Xo9XzT7pf/wVRJvd0sB8werpx\n1x4hVjYZgGpHJ41NdRQeeQ/DAMvhpOqtl/G1hQAHewPlTOhqxRvvImE4KewM8V7p+czobKLDW8An\nm3eQHw/T5sqh3eUjJxHDk4zRg4uIK4eo002Dp3fE+8r3/8TcyGFaXD72n78QC4jHXXR5vLT6JhLN\n+qBFnnPZpUQOvIOvo4XiQ7tpqJh90jS0E80q0PLNIsMhrb9Ot99+Oy0tLWzfvp1EIsG8efMIBoOZ\nrk1kVDiX67yDjWDv31LvWxjFyvETnXY+OHqvcpl5EylnF0e9fnxtTRimibezlaTThenyMDE3m67u\nXDo93RhYNATL2HjB1dQGywl2t3Fx5A1y6KE4FmFfyQxejxTiIEmrK4dAogfDglrPRM6vzGP64XY8\nTqiItZCsf48jhZV0e3LoSAQJZ+fS4i+gMnSACzxdhA8U4uv8oHu9rzXe/3KAiIyMtP5C/e///i+r\nV69m3rx5JJNJ7rvvPr73ve/xyU9+MtP1iYwprn4t9VnAwX4Dw3C6wOkk4c+l1OGgLr8YgOS7Tjzd\nYdoLJtE5sQSruJyjySQveKvoyfZTVf8egViY1px82nx5ZMe66fLmcHTiVN5esIDKhv1M6OnAFY3Q\n4/ZS7guw3+3FsJJkx7oxHS5i7mwCXe005xbyXuksDhZPB+Cmmg14mjrJbQlRP20OOR0tRIKFZwzx\ndFrj2v1MZGikFeQ//vGP+fWvf83kyb3df0eOHOH2229XkIucxqnmlfdv4U+74gr2NLR8MHUrVIdZ\nNgWSCao+ejEAh48v1+pvC5HbWg+Gk7rzqpleMQUjYeI+0NtSnhDt4I2qSzn/2C6SThctuQVUNB6k\n+sg2elzZZJs9dGflMKXpEA7TpClQRMzpJjcaJupyYxkGScNBsKuNuRV5vQXW9C3BbNEwdVbv/9QS\nFxlV0gpy0zRTIQ4wefJkkslkxooSGStOFeb9zSrOZ09TB559u1OD33rmzk8NGJsys4LDew8BBp0T\ne5eGDU2ekerezps2OdVS3mqVpFrTAPNr3sAdjzGxs4lQoAhvvIsuj49ALEwot7fFv7/kfF6Z9Qku\nf3cTF0braC8sZ+/xY9eddyH+lobeJWOP33amEFdrXGR4pRXkZWVl/Pu//zvLl/eOon3uueeYNGlS\nRgsTGSvSCfOZwRzMng7qSsohmSQ6fXbqe9GZ1RRPn41/w3/j2bef9sJyio6817sgS3eYmNdPly9I\nQ8Vsqk9jUPR4AAARhElEQVRYQrWywaLY6qLBG+BoQQXtngAToh0Ez6ugoWI2HbEoSbeHuQmTiZ0F\nNCXzsBzO1IeEhorZhCbPGNIQF5GhlVaQf+9732Pt2rU8+eSTWJbFpZdeyne/+91M1yYyZpwpzPsG\nwpXSOzVtsOlbVo6fwAVzCADxlgZym2qZWFdDtz9Iwu0h2HSUuMdLOL8kNT0sltO7f7ivJ0KwcjLR\nqguIJUwajj9+0u3pfewTlojt/2Eg3a70dENcrXGRoZVWkP/yl7/kJz/5SaZrERmVhmqZ0TOFed9A\nuGnHjzVgG9ATVkybEozha3sfCidCYzPNpRUUH9pNjy9AbkuI0OQZwAcLt3QF8mmY1rv16KmmjfWf\nGw5nbn33pxAXGTlp/XV6+eWXWbVqVWq9dZHx4lyXaB2MZZqpIDtVoPf/wNAXjn2B3re7mGffblwd\nbcRKJmH5AzjDHRTm+PHsjZDr9nA0aVFSsxNvuJ2uQD41F3x00PA2Er0r1fX/3mWVeWf1nM6mK10h\nLpIZaQV5MBjk6quvZs6cOXg8ntTtDz/8cMYKExlpgy3Req4t8xM/EKRz3bxvSdgTA71vUxLLHyCy\n8ApwZ0PCxJxYiKu5kaKuMOWHt0G0G1q8RAq9ROdcNOCxPXt3kL3zLaB3YN3ZLqd6ttfCFeIimZPW\nX6Xrrrsu03WIjDrnukTrifo+EFjJ5IAPBKcLc3PHW5hvvY4BOOcvxFU9PxWeA+aeu7MBelvprc2Y\nwXxcDgf07Md97H1ik6biam3+YJ9wSE1zc4Y7AXCF6nsH16WxrOq5DGZTiItkVtpB3tXVRXt7O5Zl\nnfkOImPEuSzReiLD5SIRDhOtOYinctqAxxqsq90yTZJ1tdDZgQUk649h9VuvfdoVVwyoaU9Dywdb\nh3a0YQaCeLK9xCZNhWwvZt7EgSHtdGEWlaamuplFJacN8Q8zEl0hLpJ5aW+a8otf/IK8vDwMw8Cy\nLAzDYMOGDZmuT2TEfdiBbpZp4vT78c6eg+F0DvrBoH+gGy4XjtIyknVHMQBHyaSTfn7AtfTifMxp\nFamlX/eVnQdOZ28LfULwpG51OH69vW+K2wkhPlRTyBTiIsMjrb9Qv/vd73jppZfIyzu7gTAicnZd\n9KlAr56Pc3Z16v5ncuLSr3z848TffgOaGzFq9+Oqnj8UTyUtCnCR4ZVWkBcVFZGbqyUZRc5VXxd9\nuvqH4ZkGxfU5cac1mnv3DrdC9UM2he50FOAiI+O07+y+fcgDgQA33HADixYtwul0pr5/++23Z7Y6\nkTGka/u2s57Kls6UtcGcuNNapkJc4S0y8k777v7d737H9ddfT3X12U1NEZGBzmUq22BT1vo7U7C7\nBtkTfSgovEVGl9O+wwOBgFrdIkPgbKeynWrKWn/pBPtQrUgnIqPXad/lWslNZOiczVS2001ZOxUF\nrsj4dNq/Dvv27eOqq6466XZNPxM5N+m2kNOZsnbiz2d6MJuIjE6nfedPnTqVf/3Xfx2uWkTkuLPp\nih/K9eBFxH5OG+RZWVnad1xkhKTTFT+U68GLiD05TvfN+fOHbxEJETnZmUK5r+UOfKj14EXEvk77\nrr/vvvuGqw4ROUdDsR68iNjXaVvkImIPCnGR8UtBLiIiYmMKchERERtTkIuIiNiYglxERMTGFOQi\nIiI2piAXERGxMQW5iIiIjSnIRUREbExBLiIiYmMKchERERtTkIuIiNiYglxERMTGFOQiImmwTHOk\nSxAZlLZMEhE5g8ibW4nX1ZJVWobv4gUjXY7IAGqRi4ichmWaxOtqAYjX1aplLqNORoN8+/btrFix\nAoD333+fG2+8kS984Qvcf//9JJPJTB5aRGRIGC4XWaVlAGSVlmnvdxl1MhbkTz31FGvWrCEajQLw\n8MMPs2rVKn79619jWRYbNmzI1KFFRIaU7+IFTLj6GnWry6iUsSCfMmUKjz32WOrrXbt2sXDhQgAW\nLVrE5s2bM3VoEZEhp5a4jFYZC/IlS5bg6nfiW5aFYRgA+Hw+Ojs7M3VoERGRcWPYBrs5HB8cKhKJ\nEAgEhuvQIiIiY9awBfns2bPZsmULAJs2bWLBAl1rEhER+bCGLci/9a1v8dhjj3HDDTcQj8dZsmTJ\ncB1aRERkzMro6I3y8nJ+85vfAFBZWcmvfvWrTB5ORERk3NGCMCI2oEVIRORUNJ9CZJTT8qAicjpq\nkYuMYloeVETOREEuMoppeVARORP9VRAZ5XwXL8AyTYW4iAxKLXIRG1CIi8ipKMhFRERsTEEuIiJi\nYwpyERERG1OQi4iI2JiCXERExMYU5CIiIjamIBcREbExBbmIiIiNKchFRERsTEEuIiJiYwpyERER\nG1OQi4iI2JiCXERExMYU5CIiIjamIBcREbExBbmIiIiNKchFRERsTEEuIiJiYwpyERERG1OQi4iI\n2JiCXERExMYU5CIiIjamIBcREbExBbmIiIiNKchFRERsTEEuIiJiYwpyERERG1OQi4iI2JiCXERE\nxMYU5CIiIjamIBcREbExBbmIiIiNKchFRERsTEEuIiJiYwpyERERG3MN9wGvu+46/H4/AOXl5Tz8\n8MPDXYKIiMiYMaxBHo1GsSyLp59+ejgPKyIiMmYNa9f63r176e7uZuXKldx8881s27ZtOA8vIiIy\n5gxrizw7O5tbbrmFz33ucxw6dIgvf/nLvPDCC7hcw97DLyIiMiYMa4JWVlYydepUDMOgsrKSYDBI\nY2MjpaWlw1mGiIjImDGsXevPPfcc3//+9wFoaGggHA5TWFg4nCWIiIiMKcPaIl++fDn33HMPN954\nI4Zh8NBDD6lbXURE5EMY1hR1u9386Ec/Gs5DioiIjGlaEEZERMTGFOQiIiI2piAXERGxMQW5iIgN\nWKY50iXIKKUh4yIio1zkza3E62rJKi3Dd/GCkS5HRhm1yEVERjHLNInX1QIQr6tVy1xOoiAXERnF\nDJeLrNIyALJKyzC09oacQGeEiMgo57t4AZZpKsRlUGqRi4jYgEJcTkVBLiIiYmMKchERERtTkIuI\niNiYglxERMTGFOQiIiI2piAXERGxMQW5iIiIjSnIRUREbExBLiIiYmMKchERERtTkIuIiNiYglxE\nRMTGFOQiIiI2piAXERGxMQW5iIiIjSnIRUREbExBLiIiYmMKchERERtTkIuIiNiYglxERMTGFOQi\nIiI2piAXERGxMQW5iIiIjSnIRUREbExBLiIiYmMKchERERtTkIuIiNiYglxERMTGFOQiIiI2piAX\nERGxMQW5iIiIjSnIRUREbExBLiIiYmOu4TxYMpnkgQce4N1338XtdvPggw8yderU4SxBRERkTBnW\nFvmf//xnYrEYzz77LN/4xjf4/ve/P5yHFxERGXOGNcjffPNNPvaxjwEwb948du7cOZyHFxERGXOG\ntWs9HA7j9/tTXzudTkzTxOU6uYxEIgFAqKF+2OoTEREZSX2Z15eB6RjWIPf7/UQikdTXyWRy0BAH\naGxsBOD2L/+fYalNRERktGhsbEx7DNmwBvn8+fN5+eWXueaaa9i2bRszZsw45c/OnTuXdevWUVhY\niNPpHMYqRURERkYikaCxsZG5c+emfR/DsiwrgzUN0Ddq/b333sOyLB566CGqqqqG6/AiIiJjzrAG\nuYiIiAwtLQgjIiJiYwpyERERG1OQi4iI2NiwjlpPl5ZyPTvXXXddan5+eXk5Dz/88AhXNLps376d\nRx99lKeffpr333+fu+++G8MwmD59Ovfffz8Ohz7PwsDXaffu3dx6661UVFQAcOONN3LNNdeMbIGj\nQDweZ/Xq1Rw7doxYLMZtt93Geeedp3PqBIO9TqWlpTqnBpFIJFizZg01NTUYhsF3vvMdPB7PWZ1T\nozLI+y/lum3bNr7//e/zL//yLyNd1qgUjUaxLIunn356pEsZlZ566inWr1+P1+sF4OGHH2bVqlVc\ncskl3HfffWzYsIHFixePcJUj78TXadeuXXzpS19i5cqVI1zZ6LJ+/XqCwSCPPPIIbW1tLFu2jJkz\nZ+qcOsFgr9NXvvIVnVODePnllwF45pln2LJlCz/+8Y+xLOuszqlR+bFRS7mmb+/evXR3d7Ny5Upu\nvvlmtm3bNtIljSpTpkzhscceS329a9cuFi5cCMCiRYvYvHnzSJU2qpz4Ou3cuZONGzdy0003sXr1\nasLh8AhWN3pcffXV3HHHHQBYloXT6dQ5NYjBXiedU4P71Kc+xdq1awGora0lEAic9Tk1KoP8VEu5\nysmys7O55ZZb+MUvfsF3vvMd7rzzTr1W/SxZsmTA6oGWZWEYBgA+n4/Ozs6RKm1UOfF1qq6u5pvf\n/Cbr1q1j8uTJPP744yNY3ejh8/nw+/2Ew2G+9rWvsWrVKp1TgxjsddI5dWoul4tvfetbrF27lmuv\nvfasz6lRGeRns5TreFdZWcnf/d3fYRgGlZWVBIPB1PK2crL+15kikQiBQGAEqxm9Fi9enFpZavHi\nxezevXuEKxo96urquPnmm1m6dCnXXnutzqlTOPF10jl1ej/4wQ948cUXuffee4lGo6nb0zmnRmWQ\nz58/n02bNgGccSnX8e65555LbQfb0NBAOBymsLBwhKsavWbPns2WLVsA2LRpEwsWLBjhikanW265\nhR07dgDw6quvMmfOnBGuaHRoampi5cqV3HXXXSxfvhzQOTWYwV4nnVODe/755/n5z38OgNfrxTAM\n5s6de1bn1Khc2U1LuaYvFotxzz33UFtbi2EY3HnnncyfP3+kyxpVjh49yte//nV+85vfUFNTw733\n3ks8HmfatGk8+OCDWsv/uP6v065du1i7di1ZWVkUFBSwdu3aAZe7xqsHH3yQP/zhD0ybNi1127e/\n/W0efPBBnVP9DPY6rVq1ikceeUTn1Am6urq45557aGpqwjRNvvzlL1NVVXVWf6dGZZCLiIhIekZl\n17qIiIikR0EuIiJiYwpyERERG1OQi4iI2JiCXERExMYU5CJjyJYtW1ixYkXaP382Pysio5OCXGQc\ne/3110e6BBH5kBTkImOcaZqsWbOGG264gauuuop/+Id/oKenhwcffBCAz33uc0DvClLLly9n2bJl\n3H777bS2tgJw5ZVX8pOf/ITly5fzmc98JrWJ0Z49e/jc5z7Htddeyxe/+EXq6+u56667ePbZZ1PH\nXrFiBdu3bx/mZywyvijIRca4t99+m6ysLJ599ln+9Kc/EY1G+ctf/sKaNWsA+M///E9aWlr40Y9+\nxC9+8Quef/55rrjiCh599NHUYwSDQZ577jn+/u//PrWc5J133sk//uM/8l//9V9cc801/Md//Aef\n/exnWb9+PQDHjh2jpaWFCy+8cPiftMg4op1IRMa4j3zkIwSDQdatW8fBgwc5dOgQXV1dA35m+/bt\nqU0uoHeZ5AkTJqS+37et8PTp0/njH/9IS0sLjY2NfPKTnwTgC1/4AtC7u9y9997L0aNH+f3vf8/S\npUuH4ymKjGsKcpExbsOGDfzsZz/j5ptv5vrrr6e1tZUTV2ZOJBLMnz+fJ598EoBoNDpgB0KPxwOQ\n2loxKytrwP2j0SihUIjJkyezbNky/ud//ocXXniBf/u3f8vkUxMR1LUuMua9+uqr/M3f/A2f/exn\nKSgo4I033iCRSADgdDoxTZMLL7yQbdu2UVNTA8ATTzzBD3/4w1M+Zm5uLiUlJbzyyisA/P73v+en\nP/0pANdffz3PPPMMJSUlFBcXZ/jZiYha5CJjzNatW7noootSX1dXV7NlyxZeeOEF3G438+bN4+jR\nowBcddVVLF26lN/97nc89NBDrFq1imQySXFxMY888shpj/PII4/wwAMP8MMf/pC8vLxU8JeWllJS\nUsJ1112XuScpIina/UxEhoxlWYRCIVasWMF///d/43a7R7okkTFPXesiMmRefPFFli5dyte//nWF\nuMgwUYtcRETExtQiFxERsTEFuYiIiI0pyEVERGxMQS4iImJjCnIREREbU5CLiIjY2P8HNyHNFPZF\nRkYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbabad30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create multi-var Gaussian model\n",
    "multi_normal = stats.multivariate_normal(mu, cov)\n",
    "\n",
    "# create a grid\n",
    "x, y = np.mgrid[0:30:0.01, 0:30:0.01]\n",
    "pos = np.dstack((x, y))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot probability density\n",
    "ax.contourf(x, y, multi_normal.pdf(pos), cmap='Blues')\n",
    "\n",
    "# plot original data points\n",
    "sns.regplot('Latency', 'Throughput',\n",
    "           data=pd.DataFrame(X, columns=['Latency', 'Throughput']), \n",
    "           fit_reg=False,\n",
    "           ax=ax,\n",
    "           scatter_kws={\"s\":10,\n",
    "                        \"alpha\":0.4})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# select threshold $\\epsilon$\n",
    "1. use training set $X$ to model the multivariate Gaussian\n",
    "2. use cross validation set $(Xval, yval)$ to find the best $\\epsilon$ by finding the best `F-score`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"../img/f1_score.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def select_threshold(X, Xval, yval):\n",
    "    \"\"\"use CV data to find the best epsilon\n",
    "    Returns:\n",
    "        e: best epsilon with the highest f-score\n",
    "        f-score: such best f-score\n",
    "    \"\"\"\n",
    "    # create multivariate model using training data\n",
    "    mu = X.mean(axis=0)\n",
    "    cov = np.cov(X.T)\n",
    "    multi_normal = stats.multivariate_normal(mu, cov)\n",
    "\n",
    "    # this is key, use CV data for fine tuning hyper parameters\n",
    "    pval = multi_normal.pdf(Xval)\n",
    "\n",
    "    # set up epsilon candidates\n",
    "    epsilon = np.linspace(np.min(pval), np.max(pval), num=10000)\n",
    "\n",
    "    # calculate f-score\n",
    "    fs = []\n",
    "    for e in epsilon:\n",
    "        y_pred = (pval <= e).astype('int')\n",
    "        fs.append(f1_score(yval, y_pred))\n",
    "\n",
    "    # find the best f-score\n",
    "    argmax_fs = np.argmax(fs)\n",
    "\n",
    "    return epsilon[argmax_fs], fs[argmax_fs]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import f1_score, classification_report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best epsilon: 3.6148577562381784e-05\n",
      "Best F-score on validation data: 1.0\n"
     ]
    }
   ],
   "source": [
    "e, fs = select_threshold(X, Xval, yval)\n",
    "print('Best epsilon: {}\\nBest F-score on validation data: {}'.format(e, fs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# visualize prediction of `Xval` using learned $\\epsilon$\n",
    "1. use CV data to find the best $\\epsilon$\n",
    "2. use all data (training + validation) to create model\n",
    "3. do the prediction on test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def select_threshold(X, Xval, yval):\n",
    "    \"\"\"use CV data to find the best epsilon\n",
    "    Returns:\n",
    "        e: best epsilon with the highest f-score\n",
    "        f-score: such best f-score\n",
    "    \"\"\"\n",
    "    # create multivariate model using training data\n",
    "    mu = X.mean(axis=0)\n",
    "    cov = np.cov(X.T)\n",
    "    multi_normal = stats.multivariate_normal(mu, cov)\n",
    "\n",
    "    # this is key, use CV data for fine tuning hyper parameters\n",
    "    pval = multi_normal.pdf(Xval)\n",
    "\n",
    "    # set up epsilon candidates\n",
    "    epsilon = np.linspace(np.min(pval), np.max(pval), num=10000)\n",
    "\n",
    "    # calculate f-score\n",
    "    fs = []\n",
    "    for e in epsilon:\n",
    "        y_pred = (pval <= e).astype('int')\n",
    "        fs.append(f1_score(yval, y_pred))\n",
    "\n",
    "    # find the best f-score\n",
    "    argmax_fs = np.argmax(fs)\n",
    "\n",
    "    return epsilon[argmax_fs], fs[argmax_fs]\n",
    "\n",
    "\n",
    "def predict(X, Xval, e, Xtest, ytest):\n",
    "    \"\"\"with optimal epsilon, combine X, Xval and predict Xtest\n",
    "    Returns:\n",
    "        multi_normal: multivariate normal model\n",
    "        y_pred: prediction of test data\n",
    "    \"\"\"\n",
    "    Xdata = np.concatenate((X, Xval), axis=0)\n",
    "\n",
    "    mu = Xdata.mean(axis=0)\n",
    "    cov = np.cov(Xdata.T)\n",
    "    multi_normal = stats.multivariate_normal(mu, cov)\n",
    "\n",
    "    # calculate probability of test data\n",
    "    pval = multi_normal.pdf(Xtest)\n",
    "    y_pred = (pval <= e).astype('int')\n",
    "\n",
    "    print(classification_report(ytest, y_pred))\n",
    "\n",
    "    return multi_normal, y_pred\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             precision    recall  f1-score   support\n",
      "\n",
      "          0       0.99      1.00      0.99       149\n",
      "          1       1.00      0.60      0.75         5\n",
      "\n",
      "avg / total       0.99      0.99      0.99       154\n",
      "\n"
     ]
    }
   ],
   "source": [
    "multi_normal, y_pred = predict(X, Xval, e, Xtest, ytest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFbCAYAAAAuggz2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPWd//H3mZnMJOQ2CbkDgYAgSgShCNi1dL2kWCyC\nFuu6Cr+HuPvw4dZqHlZUEFHBihVda1mtXdffY1ukP7TWVVq3shZRuoJBrYAJgopRhFwmIZdJQjKT\nmTm/P0LG3BmQSXKS1/MfyJk5cz4zj5l5z/dyvscwTdMUAACwJNtAFwAAAE4fQQ4AgIUR5AAAWBhB\nDgCAhRHkAABYGEEOAICFOaL1wMFgUKtWrVJpaakMw9CDDz4ol8ule+65R4ZhaOLEibr//vtls/Fb\nAgCA0xW1IN++fbskafPmzSoqKtITTzwh0zRVWFio2bNna/Xq1dq2bZsKCgqiVQIAAENe1JrDl112\nmdauXStJKisrU1JSkkpKSjRr1ixJ0ty5c7Vz585oHR4AgGEhai1ySXI4HLr77rv1xhtv6Je//KXe\neecdGYYhSYqPj1dDQ0Ov+7a0tKi4uFjp6emy2+3RLBMAgEEhGAyqqqpK+fn5io2NjWifqAa5JP38\n5z/XnXfeqR/96Efy+Xzh7U1NTUpKSup1v+LiYl1//fXRLg8AgEFn06ZNmjlzZkT3jVqQv/LKK6qs\nrNTNN9+suLg4GYah/Px8FRUVafbs2dqxY4fmzJnT6/7p6emSpF/++jfKyMyKVpkAAAwansoK3Xbz\n/wlnYCSiFuTf+973tGLFCl1//fUKBAJauXKlJkyYoPvuu0//+q//qvHjx2vevHm97t/enZ6RmaXs\nnFHRKhMAgEHnVIaUoxbkI0aM0JNPPtlt+/PPPx+tQwIAMOxwEjcAABZGkAMAYGEEOQAAFkaQAwBg\nYQQ5AAAWRpADAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAWRpADAGBhBDkAABZGkAMAYGEE\nOQAAFkaQAwBgYQQ5AAAWRpADAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAWRpADAGBhBDkA\nABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAWRpADAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAW\nRpADAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAWRpADAGBhjmg9cGtrq1auXKmjR4/K7/fr\nlltuUXZ2tm6++WaNGzdOknTddddp/vz50SoBAIAhL2pBvmXLFrndbq1fv151dXVatGiRfvzjH+vG\nG2/UsmXLonVYAACGlagF+eWXX6558+ZJkkzTlN1uV3FxsUpLS7Vt2zaNHTtWK1euVEJCQrRKAABg\nyIvaGHl8fLwSEhLU2Nio2267TYWFhZo6daruuusubdq0SWPGjNFTTz0VrcMDADAsRHWyW3l5uZYu\nXaqFCxdqwYIFKigoUH5+viSpoKBA+/fvj+bhAQAY8qIW5NXV1Vq2bJmWL1+uxYsXS5Juuukm7du3\nT5K0a9cuTZkyJVqHBwBgWIjaGPkzzzwjr9erp59+Wk8//bQk6Z577tHDDz+smJgYpaWlae3atdE6\nPAAAw0LUgnzVqlVatWpVt+2bN2+O1iEBABh2WBAGAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPI\nAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEA\nsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAw\nghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIc\nAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMEe0Hri1tVUrV67U0aNH5ff7dcstt+iss87S\nPffcI8MwNHHiRN1///2y2fgtAQDA6YpakG/ZskVut1vr169XXV2dFi1apMmTJ6uwsFCzZ8/W6tWr\ntW3bNhUUFESrBAAAhryoNYcvv/xy3X777ZIk0zRlt9tVUlKiWbNmSZLmzp2rnTt3RuvwAAAMC1EL\n8vj4eCUkJKixsVG33XabCgsLZZqmDMMI397Q0BCtwwMAMCxEdYC6vLxcS5cu1cKFC7VgwYJO4+FN\nTU1KSkqK5uEBABjyohbk1dXVWrZsmZYvX67FixdLks4991wVFRVJknbs2KGZM2dG6/AAAAwLUQvy\nZ555Rl6vV08//bSWLFmiJUuWqLCwUBs2bNC1116r1tZWzZs3L1qHBwBgWIjarPVVq1Zp1apV3bY/\n//zz0TokAADDDidxAwBgYQQ5AAAWRpADAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAWRpAD\nAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBgYQQ5AAAWRpADAGBhBDkAABZGkAMAYGEEOQAAFkaQAwBg\nYREF+aefftpt2549e854MQAA4NQ4+rrxgw8+UCgU0qpVq/Szn/1MpmlKkgKBgB544AFt3bq1X4oE\ncHpKq5r6vD0vPb6fKgEQLX0G+c6dO7V79255PB49+eSTX+/kcOjaa6+NenEA+nayoD4T+xP2wODW\nZ5D/5Cc/kSS98sorWrRoUb8UBKBn3zS0z9RxCXZgcOkzyNsVFRWpqKio2/Z169ad8YIASTJNU/X/\n/UclFcyTzekKbw/5ffK+sVXJ8xfIMIwBrLB/DFR4d2QGAjIcX39VdKyJUAcGXkRBPmvWrPD/A4GA\ntm3bpvHjx0etKKD+v/8oz5OPq2nXO8pevUY2p0shv0/la1arafe7kiT3FVcOcJXRMRjCu12o+EPJ\nUyEzI0u2/OndbifUgYEXUZBfddVVnf5evHixrrvuuqgUBEhSUsE8Ne16R02731X5mtXKXL5ClevX\nqWn3u4qfNUdJBfMGusQzaiDDu2uLu+N2eSra/vBU9Hq/du3PgUAH+ldEQd7VoUOH5PF4znQtQJjN\n6VL26jXhFvjn1yyUJMXPmhNuoQ8FZzLAd5XV9HqbEQxIkkx754984ifFctV45EvNUMOkfF2Yk9p2\nvxOhbWZktYV5RlafId4RgQ70r4g+mZMnT5ZhGOHTz1JTU3XHHXdEtTDA5nQpc/mKcIhLUubyFUMi\nxL9pgPcV2l0lflKs5I8/lGSo/pzz1TApX1JbuLuqKySbTa4ajxqDAe0qq+kW7nMmnxdxiHdUWtVE\nmAP9IKJP54EDB6JdB9BNyO9T5frOEyor16+zfIv8VEO8vXX87leebi3qk2kPa0dTgyQptrpCjRMm\ny7Q7lHDogGKrKiSZqj9nuky7o+3+NW29be3h/q7HG3689hZ7pGidA9EX0beC1+vVhg0b9O6778rh\ncGju3Lm65ZZbFBsbG+36MEx1nNgWP2tOpzHy8jWrLRnmp9MKL35rh1w1Htmbj2tk3IhwK7k3RjDQ\nKexNu0O+tCzFVpVLMtSSlhW+n6vGo5bMHBmhkBonTA7f3+8eKWfdMflSM7r9cGjvCTidQCfMgeiI\nKMiXL1+u8ePH67HHHpNpmvrDH/6ge++9V48//ni068Mw5X1jazjE20O745i5942tlpq1fqohvqus\nRkYwoJE1HikUUlz5YTWNmxRuJffUMk/8pFiHPy5VfVKaDo86W7ZQUJIUsuXIdnambKGgco6WKvnj\nP6s+KU2lilOyt1r1SWkad+LxEj8pbgtx90g1TMoPj61LncfXTyfQCXMgOiIK8qNHj+rXv/51+O97\n771XP/jBD6JWFJA8f4EkdTqPvD3M288jt4pTCfGOY9+m3SFfaoZcNR41Z+fKtNm6tZJ3fNJ2f1so\nqCkHSyVJyd5qjQ0FNe7Ix5IMleZOkSS56zxK8Vapxp2pZG+1Ss6eo7LMPAUcTh3+pKbTY0yUpAN7\nlXxwn5x1x+R3p3UaX+9ab6SBTpgDZ15EQT527Fi9//77mjlzpqS2MfOxY8dGtTAMb4Zh9Njitjld\nQ7Yl3tMEtoZJ+eEWuOFvkemMDYd3RyGbXfVJaUr2VsubkKJk7zGNaG6UaRhKqauUKUOmzSbJlGGG\nVJecodHln8ld71FdcoYOjzq702N80ORUcvkXClZWaYy/TqYjptP4ek+1E+bAwIgoyA8fPqwbbrhB\neXl5stvtKi0tVXJysi655BIZhqFt27ZFu07AciIN8ZPNQDftDn2xfVe4G1yjzu7xfodHnS1b9lkK\n2eyaufcvGtHcIJ8zVnVJaTJtdiV7q1Wam68j2WdJkr5T9KriWhqVUl+lIyf2OzzqbDnTc+V3xin3\n6EGleD2qCLTKG3Ap4cT4el/PgzAH+l9EQf7MM89Euw5gSDlTIS5J/3ugSlO81ZLaus3bw7orWyio\nfSWVsoeCSqlu1RfOccpsPqbGg1+oIi5VB5PyFKyzS3XVcgb8mlzdpMRAs2Kz4mULBcM/ANJqylSd\nmqP3p10mIxRUsveY6pNGyjzaoOSP/6zcc/J6nXBHmAP9L6Igz8jI0M6dO1VbW9tpOxdSAbr7piHe\nrQu9Q5d3fVJaOMTbw7f4ozKd7T2sabWf6e989TrmTJIr1CpX0K8YM6gyI0MZLbX6LGGUJOls72Fl\ntNSqzpmoirhUxXpblfj2NtXFxMt+7JDkjlNaTZmc/mYlNbZ95pMbjklty0jo8MelSu2li739eZ3q\nrHYApy+iIL/99ttVVVWlCRMmdLpQBUEOdPZNQzxz+580ouxLlTjTpWmXhbd37DaXpNZ33lVGS63i\nAj6l2mOU3XxMia3NSvF7FZBdpiF9mjhKmS11spmmXEG/fnj4LRkyZUg6OiJdzQ6XilLP0eyaj5V9\nvFqTAofVarOrtq5ZX43IVMDhlMvfrLSaMh1zZ8m02cKt9bcPeTV3Umq30906Pr9IwpxWOfDNRRTk\nn3/+uV5//fVo1wJY2ploiSd98pGa6hs1NuaY9kyZq4DDGb49ZLNr74mu8++01MpmmsptqtSniTmS\naaoxJk4xZpLqnfGqcyYqZHNob+pEfZYwSn9XtU+JgeNKbm2SYZoyJe1JnaRmZ5yqXcma5D2shph4\nlceNVFHqOco7XqGsN7cqPr5ZLa4RSqst1/HYESrLzJNp2GQLBfXF9l2aGdfc67nttMyB/hFRkOfm\n5qqsrEw5OTnRrgcY0vocE7c7VBmwK0GSzxnbrQt9b0mlJClos6valazz6j7XSF+dcpqrVBWbor+m\nT1VxyoTwfapKShQ0WiQd0sjWJiW3tigu2CRP+rgTY+a5kqTilAkyJI301csTmyK/w6mMllqFDEPH\nqps0Oj6k43GJkgzJlOqT0yS1jdd/6pUmqvdz2yNBqxz4Zvr85C1ZskSGYaimpkYLFizQ5MmTZbd/\nPcnmt7/9bdQLBKzgTFz85O1DXuWe+x2l1FSoNjVLIZtduUcPKtlbrb/V2CV3niTpvNpDSm+pky0U\nUKvNoRR/ozJaajXt812qjjmqPTEZCho2ybCFH3ufI037HGk6v9WjnLomHWlo0tHaT2Q3Q0qfMkUf\npUyQ/cQCMkGbXZ7YFGW01GpfylnaK+l8t6n6pLTw7HZJ4XH7jue2d+1qp1UORF+fQf6Tn/ykv+oA\nhry+WuPtE9sOjzo7HJa2UFDJ3modKfMqQ9LBUK7Oqf9CF1Z9pIaY+LagNk2FZCjheL1i7Sm6oLVS\nWaHjKrMn6MOYDEnSt1orlBFqVoUtXuNDXuUF6hSrttDOCjWpYm+1POddFJ4E54lN0X53nmy1IaWd\naKX/v0Cu8kd17pFrH7f/SHbNVfcrqQHoH30G+axZsyRJ7733XqfthmHI5XLJ6/UqKSkpetUBw1B7\nizdks+tvNXZlSPLEpsgeCirNV68GxwgltjZpV/p52uF16OJQSBmSRocalGE2K2RIIcOQ3Qzp/FaP\nLmn9Sl7DJYcZ1PjWWoVsNuUF6tQimwKGXRe0lqt8z591VopL5XEj22a4B0YpzVcvScpsPqaDSbnh\nLv72fzvW2tPFVk63qx3AqYnok/bUU0+puLhYF154oUzT1O7duzVq1Cg1Njbq9ttv73W51r179+qx\nxx7Txo0btX//ft18880aN26cJOm6667T/Pnzz9gTAQZKJN3qkbTGe7LfnaeDoVyd7T2sb1cXKy7g\nU/mINO13JuqjlAmqKPdrQqBeObZG+WVXg+GUO+TTB444SVKG2ax6w6lk06daORVjmIoPtehvMZk6\nYk9STrBRhgwFDJuajh2TbVSqKuJS5Xc45YlN0bSaTyXD0PfKdqvZ49KEtBj5nHHh9dw7al9OtqeL\nrQCInog+baZpasuWLeHJbpWVlVq5cqU2btyoJUuW9Bjkzz77rLZs2aK4uLYvlJKSEt14441atmzZ\nGSwfGDo6tnQlhSe32UNBZbS0nc/dYo9RrTNRqb56Ze99W1WONNXbY3VAqRobbNCX9kS9F5OlD2La\nVmHz2OJkM019aMQqTT59GJMhRyioP7kmKGjYtM+RpqmBamWFmvReTLY8mecreKKGg0m5ymqukSlT\nExuO6FDCKKXVVKksIy+8MI0kjS7/TCNr2mavH7tgLiEO9LOIPnEej6fTjPXMzEx5PB4lJCTINM0e\n98nNzdWGDRt01113SZKKi4tVWlqqbdu2aezYsVq5cqUSEhLOwFMArK99Ulv7hLJ259aVKrP5mGKD\nrXL7G2TI1JS6UvntMUpsbdE+R5oqbG0zvt+LydIex4mJbpKmt3r0rdZK2WSqLCZBCcFGzWitVJMR\no8OB5PAY+ocxGbKbIQUNm7JsdtlDQQVtdgVtdlXEpSqjpVblsSPbZsun5si02VSflBZeqz2l3iOd\nPUGuE2uxA+hfEQX59OnT9dOf/lQLFixQKBTSa6+9punTp+utt97SiBEjetxn3rx5OnLkSPjvqVOn\n6pprrlF+fr5+9atf6amnntLdd999Zp4FYGHtk9okKe9wifK+3CcZdgWCWcqvPaTEwHE12VyqGJEm\n07BpVHO1Wm0OGWpbnOnDmAw5QwH5bV9/nJ2hgHKCjUo2/TJMU7mBerXIrlrDJdMwlBNqlM0fUobZ\nrGojVu85syW1/XDoOOFtvztP474qU3bLMQUNm96f9v3w5VGnHHxXpmGTZCi2/CvJZlfCoQNMdAP6\nWURBvmbNGm3evFkvvPCC7Ha7vv3tb+tHP/qR3nnnHT366KMRHaigoCA8Ma6goEBr1649/aqBIaT9\nqmPueo+MUEix/mZJUkagRoZMJfsblWw0qS42SU2OOBUnj5ff5tD+luMKGjZNb/UoJ9SoMlvbTPX2\nWeojFJDLDCgrdFxOhfSJw606m0uGDFXYRigj1KxvtVYqzWxWXtCrP8RODHfhZ7TU6mAoV/ZQUKOb\nqyRJo5urVBnwhxepaT/97IvRk5UW1yLTZmOiGzAAIu5av+SSS3TJJZd02vbd73434gPddNNNuu++\n+zR16lTt2rVLU6ZMOfVqgSGq/bSz0eWfKe9wsSRD5aE0VcalKrG1WV5nvI7bXdqZlq+zGo8qs/mY\npLaW96zWCiWbLRptNMoRCujvA0flNVwqs43QR/Y0+Y1apYR8ygvU682YMfqbM0tBw6YL/OVKM5vV\nLIfGhrxyhgLh88c9sSnh7vWvRmRqzPFKBcdP6LTSXMdlY/NVxkQ3YIBE9Im74YYbwmust7a2qrq6\nWuecc47+8Ic/RHygBx54QGvXrlVMTIzS0tJokQNdtF9GtH2MPObEam6GpPzaQ8pqqdXZDV8praVW\npmHT3yf6NLruS00M1MphhOQyAsqwJajecCovWC/DNHXMHiev4VSu2SBJmhGs1t/UNhHuPWe28oJe\njQ155ZddC9OD8kj6a8a08IQ3Sfrz6AvlDPh1zrQxPdY8d1KqGpTaY0ucxWCA6IsoyN98881Of+/b\nt0+bNm066X6jR4/Wiy++KEmaMmWKNm/efBolAsNL10uU7k8ep/SWOoUMQ/l1n8swQ5Ik07ApITtL\nrYfLVGXEqs5wqcI+QjaZsoVC+sqeqJBh6E+uc7XU3K8ktcrU15NT7WZIL8adrbigXwvTTywQ01wT\nXrq1o55CvKvTbYmzPCvwzdhOfpfupk6dqpKSkjNdCzBk9dUynTup99umTckMzx63maZkmiobka6K\nuDR95B4v07Dpk1Hnqd6IlQybJgTrNCrUqLNMry5qPaoEs1XNdqd2O3O03z5S78VkK2jY9C1/hb7n\n/1KX5MQqdUq+PLEpym4+pszmYzrbe7hbDadTO61xoH9E9BP63/7t3zr9/dlnn2nkyJFRKQiwmrz0\n+NNea723y4B2NP2cNO0rCepgUm54GdWKuFTtd+eF10f/rt2p7OYqTW04Kl+MU7bjIVUbcTouh+xm\nSB/GZGifI03pU6ZoXs2nml3dtkKcUfOJso9XyxOXqorYFJmGLTzRLWizh0O86znuUt8hHila48A3\nd1p9YRdccIGuuOKKM10LMKRdmJPaaYW3xAN75ao7Jl9qhuZOyu9xhbfcoweVd7hY35Ghv4Rywiu9\ntY9ht/8bF2zR+MYyBQ27TMNQU3K6fM5kNadOVPqJwA/a7Dqv9pBmV5fIGQpodFOlEoItGtVcrSO+\nDH3kHh9eW71jiHc8x/3wqLN7DPWeniuA/hFRkN96662qqanR3r17FQwGdf7558vtdke7NmDISjqw\nV6l/e0eB+ERJUmMwoLmTUjuFuS0UlLvOo7iWttb+tzJb5Zqcpg8/ru70WM6AX02OOB1KGKXE1uPa\nlZ6vEvd4SW1B335ueLUrWSN99WqIiVeSv0ne2FTJV6vYoE82M6SPk8dJkvLPy9G0DjW0n+Oe7K3W\n2FBQ34r3t10YRT2HdaQhTmscODMiGiP/61//qoULF+rll1/Wf/3Xf+nKK6/U9u3bo10bYBmRhtKF\nOakyggE5644pEJ8kR1ODfO6R4e71jt3VIZtdde4MNcfGqzk2QXXJGQqdaCm3t5bPrSvV31XtU1zA\nJ1OGGmNGKGBzhE8d67i8a5qvXsdcySqPG6kj8elqdTjVbHfpUMIofZg6Sfnn5Sj/vM5XOGs/x12S\nvAkp+la8X1LbhVGMYKDH53cmXy8AJxdRi/yJJ57Q7373O40Z0zZz9auvvtKtt96qiy++OKrFAUPR\nnDEZKj7Utjyq332eGiZP63T73EltYf/2IW+n09G6dmdPPydN3yn6X8XZm/RVs1QZmxoe4/4sMEp+\nhzN8bfGs5prwuLoz4Ne3q4uVkpupejOkg5NmK6bD+eFdtZ8vftHkdPk6XKqUU82AwSGiIA8EAuEQ\nl6QxY8YoFApFrSjAinqa9GYGAjIc3T9mUy76torKa3qc6Jb4SbFiqyt0RVqWXlPOScaj29Z3SEtL\n0IhYp0bWVShgj9Fo5+Hwuu2jy+vkrmtVvDtFMaPaWvKJR49L3mrVJWd0WuSlJx17CRom5X/j88Vp\njQNnVkRBnpOTo//8z//U4sWLJUkvvfSSRo0aFdXCACsal+LSF7U+SVKo+EPJUyEzI0u2/Onh+7Rv\nn5WRpaLUsZ32N4IBJX+8R44mr1xVFfruDybLtDt6nAgXstlVmjtF7nqP6hNHKqmxVuXpY5Vd9aWM\nUEh5h0uUUlMhd+Mx1bgzw1csa194xpGZ12uI9zUjnRAHBpeIgvxnP/uZ1q5dq2eeeUamaWrOnDla\ns2ZNtGsDLKXpb++rtbxMGdk5qsyeKHkq2m7wVIRb5mYg0Gn7nMnn6V2Pt8tpaGaXfzsHa8dQb+96\nD9ns4dnl1ak5JxrqpkIOhyRThhkKj7FL3Weidz1GR72dIkeAA4NDREH+29/+Vr/4xS+iXQtgWWYg\noNbyMklSa3mZ8qaer88zstpCOyMr3L1uOBwyu2yfXfOl5KnQfsWqYfI01Z8zXa7qCvnSsnoM0K6B\n2x7sHdc+t4WCGl3+mZK91SrNzQ+HvfT1TPSJmfGSmpU3IanXc9kTO4yJt1/V7FTHwglxILoiCvLt\n27ersLAwvN46gM4Mh0Mx2TlqLS9TTHaODIdDEy6+SJ+X13cbI7flT+/WQjePHtbkBq+MkYnaPSlf\njRMmR7zkaa/d4JPTZQQDSrE7NK7LTYm2POkkFzkxggG5ajyS2mapT8lI6nG8vy+EOBB9EX0q3W63\nLr/8ck2ZMkUulyu8fd26dVErDLCa+Bkzu01uG5+d3OOqbx1b6KGR6TIPFEuJSTKOVWnOlPPDt3dc\nQOZ09BbSvU1a67rvhEkTuvUqRIIAB/pPRJ/Mq666Ktp1AENCT2HXHmq9LeNqnzZTQUnGsapugdmx\nG/ubhnpXvYV4p67znNReZ973hhAH+lfEQX78+HHV19fLNM2T7wCgm77WZLdP696a76pjwJqBgN71\neL9xTZGMd0ca4gQ4MDAivmjKc889p5SUFBmGIdM0ZRiGtm3bFu36gCGlr9Z5pIHZfvra7C6ntQ0U\nAhwYWBF9c7z88st68803lZKSEu16gGGhY/idypXTup6+dqrd3mcK4Q0MHhF9A2RkZCgxMTHatQDD\nRscAPtkYekc9nb7WXwhvYHDq81ug/TrkSUlJuvbaazV37lzZ7V8vF3nrrbdGtzpgCGpfOCYmO0fx\nM2aGt0faSu94+lq0Ed7A4NfnN8HLL7+sq6++WlOnTu2veoAhrevCMb0Fck8BWlrVFL5/NEKc0Aas\nqc9vg6SkJFrdwBnU08Ixkcr46uMeW/Id9daSJ6SBoavPbxFWcgPOvJ4WjjmZb9KSBzC09flN8umn\nn+rSSy/ttp3Tz4Bv5lS7xr9JSx7A0Nbnt8HYsWP17//+7/1VC4A+nE5LHsDQ1+c3QkxMDNcdBwYR\nQhxAV7a+bpwxY0Z/1QGgB2YgMNAlABjk+vx5v3r16v6qA0AXvZ1vDgAd9dkiBzAwepqlDgA9IciB\nQah9lrokZqkD6BPfDsAgxSx1AJGgRQ4MYoQ4gJMhyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBQG0X\ng6p7bYtCfl+n7SG/T3WvbZFpmgNUGdA3ghwAJNX/9x/lefJxla9ZHQ7zkN+n8jWr5XnycdX/9x8H\nuEKgZwS9tLToAAANMElEQVQ5AEhKKpin+Flz1LT7XZWvWa1AfZ3K16xW0+53FT9rjpIK5g10iUCP\nCHIAkGRzupS9ek04zD+/ZmE4xLNXr5HN6RroEoEeEeQAcILN6VLm8hWdtmUuX0GIY1AjyAHghJDf\np8r16zptq1y/rtsEOGAwIcgBQF9PbGvvTh//+1c7jZkT5hisCHIAkOR9Y2unMXFHsrvTmLn3ja0D\nXSLQo6hekWHv3r167LHHtHHjRn355Ze65557ZBiGJk6cqPvvv182G78jAAwOyfMXSGqbvd4+Jt4+\nAc77xtbw7cBgE7UkffbZZ7Vq1Sr5fG3dUevWrVNhYaF+97vfyTRNbdu2LVqHBoBTZhiG3Fdc2W1i\nm83pkvuKK2UYxgBVBvQtakGem5urDRs2hP8uKSnRrFmzJElz587Vzp07o3VoAACGjagF+bx58+To\ncC1l0zTDv2jj4+PV0NAQrUMDADBs9Nsgdcfx8KamJiUlJfXXoQEAGLL6LcjPPfdcFRUVSZJ27Nih\nmTNn9tehAQAYsvotyO+++25t2LBB1157rVpbWzVvHusWAwDwTUX19LPRo0frxRdflCTl5eXp+eef\nj+bhAAAYdjiRGwAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAAL\nI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPI\nAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEA\nsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACxv0QW6apupe26KQ39dpe8jv\nU91rW2Sa5gBVBgDAwHMMdAEn0/DWNjX+5v+qadc7yl69RjanSyG/T+VrVqtp97uSJPcVVw5wlQAA\nDIxB3yJPvGiu4mfNUdPud1W+ZrUC9XXhEI+fNUdJBfMGukQAAAbMoA9yI8ap7NVrwmH++TULwyHe\n3kIHAGC46veu9auuukoJCQmSpNGjR2vdunUn3cfmdClz+Qp9fs3C8LbM5SsIcQDAsNevQe7z+WSa\npjZu3HhK+4X8PlWu7xz4levX0SIHAAx7/dq1fuDAATU3N2vZsmVaunSp9uzZc9J9zFZ/pzHx8b9/\ntdOYedfZ7AAADCf9GuSxsbG66aab9Nxzz+nBBx/UnXfeqUAg0Oc+Df+7o9OYuCPZ3WnM3PvG1n6q\nHgCAwadfu9bz8vI0duxYGYahvLw8ud1uVVVVKTs7u9d9Ev/+Uo1wpyipYF64G93mdCl79Rp539iq\n5PkL+qt8AAAGnX5tkb/00kt65JFHJEmVlZVqbGxUenp6n/sYhiH3FVd2Gwu3OV1yX3GlDMOIWr0A\nAAx2/doiX7x4sVasWKHrrrtOhmHo4YcflsMx6NekAQBg0OrXFHU6nXr88cf785AAAAxpg35BGAAA\n0DuCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAw\nghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIc\nAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAA\nCyPIAQCwMIIcAAALI8gBALAwghwAAAsjyAEAsDCCHAAACyPIAQCwMIIcAAALI8gBALAwghwAAAtz\n9OfBQqGQHnjgAR08eFBOp1MPPfSQxo4d258lAAAwpPRri/wvf/mL/H6/XnjhBf30pz/VI4880p+H\nBwBgyOnXIP/ggw/0ne98R5J0/vnnq7i4uD8PDwDAkNOvXeuNjY1KSEgI/2232xUIBORwdC8jGAxK\nkjyVFf1WHwAAA6k989ozMBL9GuQJCQlqamoK/x0KhXoMcUmqqqqSJN128//pl9oAABgsqqqqIp5D\n1q9BPmPGDG3fvl3z58/Xnj17NGnSpF7vm5+fr02bNik9PV12u70fqwQAYGAEg0FVVVUpPz8/4n0M\n0zTNKNbUSfus9U8++USmaerhhx/WhAkT+uvwAAAMOf0a5AAA4MxiQRgAACyMIAcAwMIIcgAALKxf\nZ61HiqVcT81VV10VPj9/9OjRWrdu3QBXNLjs3btXjz32mDZu3Kgvv/xS99xzjwzD0MSJE3X//ffL\nZuP3rNT5ddq/f79uvvlmjRs3TpJ03XXXaf78+QNb4CDQ2tqqlStX6ujRo/L7/brlllt01lln8Z7q\noqfXKTs7m/dUD4LBoFatWqXS0lIZhqEHH3xQLpfrlN5TgzLIOy7lumfPHj3yyCP61a9+NdBlDUo+\nn0+maWrjxo0DXcqg9Oyzz2rLli2Ki4uTJK1bt06FhYWaPXu2Vq9erW3btqmgoGCAqxx4XV+nkpIS\n3XjjjVq2bNkAVza4bNmyRW63W+vXr1ddXZ0WLVqkyZMn857qoqfX6cc//jHvqR5s375dkrR582YV\nFRXpiSeekGmap/SeGpQ/G1nKNXIHDhxQc3Ozli1bpqVLl2rPnj0DXdKgkpubqw0bNoT/Likp0axZ\nsyRJc+fO1c6dOweqtEGl6+tUXFyst956S9dff71WrlypxsbGAaxu8Lj88st1++23S5JM05Tdbuc9\n1YOeXifeUz277LLLtHbtWklSWVmZkpKSTvk9NSiDvLelXNFdbGysbrrpJj333HN68MEHdeedd/Ja\ndTBv3rxOqweapinDMCRJ8fHxamhoGKjSBpWur9PUqVN11113adOmTRozZoyeeuqpAaxu8IiPj1dC\nQoIaGxt12223qbCwkPdUD3p6nXhP9c7hcOjuu+/W2rVrtWDBglN+Tw3KID+VpVyHu7y8PF155ZUy\nDEN5eXlyu93h5W3RXcdxpqamJiUlJQ1gNYNXQUFBeGWpgoIC7d+/f4ArGjzKy8u1dOlSLVy4UAsW\nLOA91YuurxPvqb79/Oc/19atW3XffffJ5/OFt0fynhqUQT5jxgzt2LFDkk66lOtw99JLL4UvB1tZ\nWanGxkalp6cPcFWD17nnnquioiJJ0o4dOzRz5swBrmhwuummm7Rv3z5J0q5duzRlypQBrmhwqK6u\n1rJly7R8+XItXrxYEu+pnvT0OvGe6tkrr7yiX//615KkuLg4GYah/Pz8U3pPDcqV3VjKNXJ+v18r\nVqxQWVmZDMPQnXfeqRkzZgx0WYPKkSNHdMcdd+jFF19UaWmp7rvvPrW2tmr8+PF66KGHWMv/hI6v\nU0lJidauXauYmBilpaVp7dq1nYa7hquHHnpIf/7znzV+/PjwtnvvvVcPPfQQ76kOenqdCgsLtX79\net5TXRw/flwrVqxQdXW1AoGA/vmf/1kTJkw4pe+pQRnkAAAgMoOyax0AAESGIAcAwMIIcgAALIwg\nBwDAwghyAAAsjCAHhpCioiItWbIk4vufyn0BDE4EOTCM7d69e6BLAPANEeTAEBcIBLRq1Spde+21\nuvTSS/VP//RPamlp0UMPPSRJuuaaayS1rSC1ePFiLVq0SLfeeqtqa2slSZdccol+8YtfaPHixbri\niivCFzH6+OOPdc0112jBggW64YYbVFFRoeXLl+uFF14IH3vJkiXau3dvPz9jYHghyIEh7sMPP1RM\nTIxeeOEFvfHGG/L5fHr77be1atUqSdLvf/971dTU6PHHH9dzzz2nV155RRdddJEee+yx8GO43W69\n9NJL+od/+IfwcpJ33nmn/uVf/kV//OMfNX/+fP3mN7/RD3/4Q23ZskWSdPToUdXU1GjatGn9/6SB\nYYQrkQBD3AUXXCC3261Nmzbp888/1xdffKHjx493us/evXvDF7mQ2pZJTk5ODt/eflnhiRMn6n/+\n539UU1OjqqoqXXzxxZKkf/zHf5TUdnW5++67T0eOHNGrr76qhQsX9sdTBIY1ghwY4rZt26Zf/vKX\nWrp0qa6++mrV1taq68rMwWBQM2bM0DPPPCNJ8vl8na5A6HK5JCl8acWYmJhO+/t8Pnk8Ho0ZM0aL\nFi3Sa6+9ptdff13/8R//Ec2nBkB0rQND3q5du/T9739fP/zhD5WWlqb33ntPwWBQkmS32xUIBDRt\n2jTt2bNHpaWlkqSnn35ajz76aK+PmZiYqKysLL3zzjuSpFdffVVPPvmkJOnqq6/W5s2blZWVpczM\nzCg/OwC0yIEh5v3339f06dPDf0+dOlVFRUV6/fXX5XQ6df755+vIkSOSpEsvvVQLFy7Uyy+/rIcf\nfliFhYUKhULKzMzU+vXr+zzO+vXr9cADD+jRRx9VSkpKOPizs7OVlZWlq666KnpPEkAYVz8DcMaY\npimPx6MlS5boT3/6k5xO50CXBAx5dK0DOGO2bt2qhQsX6o477iDEgX5CixwAAAujRQ4AgIUR5AAA\nWBhBDgCAhRHkAABYGEEOAICFEeQAAFjY/wdc7nH9+CU6aAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa9f7588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# construct test DataFrame\n",
    "data = pd.DataFrame(Xtest, columns=['Latency', 'Throughput'])\n",
    "data['y_pred'] = y_pred\n",
    "\n",
    "# create a grid for graphing\n",
    "x, y = np.mgrid[0:30:0.01, 0:30:0.01]\n",
    "pos = np.dstack((x, y))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot probability density\n",
    "ax.contourf(x, y, multi_normal.pdf(pos), cmap='Blues')\n",
    "\n",
    "# plot original Xval points\n",
    "sns.regplot('Latency', 'Throughput',\n",
    "            data=data,\n",
    "            fit_reg=False,\n",
    "            ax=ax,\n",
    "            scatter_kws={\"s\":10,\n",
    "                         \"alpha\":0.4})\n",
    "\n",
    "# mark the predicted anamoly of CV data. We should have a test set for this...\n",
    "anamoly_data = data[data['y_pred']==1]\n",
    "ax.scatter(anamoly_data['Latency'], anamoly_data['Throughput'], marker='x', s=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# high dimension data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mat = sio.loadmat('./data/ex8data2.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = mat.get('X')\n",
    "Xval, Xtest, yval, ytest = train_test_split(mat.get('Xval'),\n",
    "                                            mat.get('yval').ravel(),\n",
    "                                            test_size=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best epsilon: 5.240023323411777e-19\n",
      "Best F-score on validation data: 0.8571428571428571\n"
     ]
    }
   ],
   "source": [
    "e, fs = select_threshold(X, Xval, yval)\n",
    "print('Best epsilon: {}\\nBest F-score on validation data: {}'.format(e, fs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             precision    recall  f1-score   support\n",
      "\n",
      "          0       0.95      0.91      0.93        43\n",
      "          1       0.56      0.71      0.63         7\n",
      "\n",
      "avg / total       0.90      0.88      0.89        50\n",
      "\n"
     ]
    }
   ],
   "source": [
    "multi_normal, y_pred = predict(X, Xval, e, Xtest, ytest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "find 9 anamolies\n"
     ]
    }
   ],
   "source": [
    "print('find {} anamolies'.format(y_pred.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The huge difference between my result, and the official `117` anamolies in the ex8 is due to:\n",
    "1. my use of **multivariate Gaussian**\n",
    "2. I split data very differently"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
